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FOREWORD 


Nonconvex  programming  computer  programs  are  an  essential 
part  of  the  effective  practice  of  operations  research  as  applied 
to  military,  industrial,  and  economic  problems.  Many  such  pro- 
grams, however,  fail  to  converge,  find  only  local  optima,  or 
become  unstable  when  applied  to  large  problems. 

This  paper  documents  a computer  program  that  can  be  applied 
to  a broad  range  of  nonconvex  programming  problems.  The  program 
is  important  in  that  it  finds  a global  optimum  in  a finite  number 
of  steps,  and  has  proven  to  be  stable  for  large  problems. 
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A.  INTRODUCTION 

Mathematical  programming,  a fundamental  tool  of  operations 
research,  is  frequently  used  to  find  solutions  to  optimization 
problems  arising  in  the  analysis  of  military,  industrial  and 
economic  models.  The  utility  of  linear  programming,  applicable 
to  models  in  which  all  equations  are  linear,  is  well  known  and 
one  reason  for  the  widespread  use  of  linear  programming  is  the 
availability  of  computer  codes  for  solving  linear  programming 
problems.  Many  important  problems,  however,  cannot  be  conven- 
iently modelled  in  a linear  framework.  A brief  survey  of  recent 
literature  reveals  nonlinear  programming  applications  to  missile 
allocation,  failure  diagnosis,  media  selection  for  advertising, 
facility  location,  chemical  process  scheduling,  design  of 
sewers,  and  so  forth.  For  these  types  of  analyses,  nonlinear 
optimization  problems  must  be  solved. 

A major  difficulty  that  arises  in  nonlinear  programming  is 
the  existence  of  local  optima.  Except  when  certain  convexity 
conditions  obtain,  nonlinear  programming  codes  in  general  cannot 
guarantee  that  the  answers  they  produce  are  globally  optimal. 

► Although  the  use  of  local  optima  may  be  useful  in  some  cases, 

basing  analyses  on  local  optima  rather  than  global  optima 
defeats  the  purpose  of  engaging  in  mathematical  programming. 

It  is  therefore  noteworthy  when  a computer  code  becomes 
available  that  can  guarantee  a global  optimum  for  a large  class 
of  nonlinear  programming  problems — the  class  of  separable, 
piecewise  linear  problems — in  a finite  number  of  steps.  Further, 
the  code  can  generate  piecewise  linear  approximations  to  any 
* separable,  continuous  optimization  problem  and  find  a globally 

optimal  solution  of  the  approximate  problem.  The  code  has  been 
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tested  on  a wide  range  of  problems,  and  the  size  of  problems 
that  can  be  handled  is  limited  only  by  computer  storage  and  run 
time  considerations. 

The  code  is  based  on  an  algorithm  by  James  E.  Falk  of  The 
George  Washington  University.  A theoretical  treatment  of  this 
algorithm  is  reprinted  in  Appendix  A,  which  also  describes  some 
of  the  background  of  this  approach.  The  algorithm  uses  branch- 
and-bound  to  generate  a sequence  of  linear  programming  sub- 
problems. 

An  earlier  realization  of  this  algorithm,  the  NUGLOBAL  code1, 
was  found  to  have  serious  stability  deficiencies  in  its  linear 
programming  subsection  when  applied  to  large  problems.  The 
code  described  in  this  paper,  which  is  embodied  in  a program 
named  MOGG,  was  therefore  developed  to  be  stable  and  also  to 
correct  some  other,  less  serious,  computational  inefficiencies. 

In  particular,  a linear  programming  package  designed  by  John  A. 
Tomlin  of  Stanford  and  adapted  by  Paul  F.  McCoy  of  IDA  was 
incorporated  into  the  new  code.  This  linear  programming 
package  has  proved  to  be  trustworthy. 

This  paper  is  divided  into  seven  parts:  Part  B concisely 
describes  the  types  of  problems  the  code  will  solve,  and  the 
details  of  computing  piecewise  linear  approximations  to  separable, 
continuous  optimization  problems;  Part  C is  a user’s  guide 
explaining  the  input  necessary  to  run  Program  MOGG;  a sample 
problem  appears  in  Part  D;  Part  E presents  a flowchart  of  the 
algorithm  as  implemented  in  this  code;  Part  F remarks  on  some 
of  the  error  messages  that  may  be  encountered  during  a MOGG  run; 
and  Part  G comments  on  some  of  the  important  variables  and 
tolerances  used  by  the  code. 


Coffman,  Karla  R,  NUGLOBAL--User ' s Guide , Technical  Memorandum 
Serial  TM-64866,  The  George  Washington  University  Program  in 
Logistics,  Washington,  D.C.,  March  1975* 
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Appendix  A has  already  been  described.  Appendix  B is  a 
description  of  the  linear  programming  package  and  Appendix  C 
contains  a complete  FORTRAN  listing  of  Program  MOGG. 


B.  PROBLEMS  TO  WHICH  MOGG  APPLIES 

Consider  the  problem  P-1 


C minimize  F^Cx) 


where 


x = (x. 


P-1  < 


subject  to  F^(x)  < b^ 


F1(x)  = bi 


l . < x . < u . 
J ~ j " J 


i=2, . . . ,q 
i=q+l , . . . ,m 
j — 1 > • • • j ft  ■ 


We  will  assume  that  all  F^(x)  are  continuous  over  the 


rectangle  < x.  < u^  j=l,...,n  (this  condition  is  actually 


somewhat  stronger  than  necessary,  see  Appendix  A).  With  no 
further  restrictions,  this  problem  in  general  cannot  be 


solved.  However,  if  each  F^Cx)  is  separable,  i.e.,  if  each 


F^Cx)  can  be  written 


F (x)  = l F (x  ), 
j = l J J 


then  we  can  approximate  Problem  P-1  by  a piecewise  linear 
problem  in  the  following  manner.  Consider  Figure  1 which,  we 


will  imagine,  depicts  some  F..(x.)  for  l . < x.  < u..  Let  us 

^-J  J J J J 


divide  the  interval  u.]  into  t intervals  by  specifying  the 

J J 


0 1 4-  V ^ 

points  tZj,  Zj  , . . . , z j } j which  we  shall  call  "cuts"  where  all 
that  we  require  is 


i . = z.  < z.  <•••<  z 


t.  = 


j 


V 


Now  we  define  a new  function  F. .(z.)  for  i . < z.  < u.  as 

— J J J J J 

follows : 


I 

t 


Figure  1 . A TYPICAL  F.  .(x  .) 

■ 3 J 


W = zkii_^k  ^Fij(zj  } ‘ Fij(zj}j  + Fij(zj} 

for  z.  e [Zj,  zj+1]»  k=0, . . . , t-1. 

It  is  easy  to  see  that  F^.tz. ) is  continuous  and  piecewise 
linear.  Figure  2 shows  the  approximation  F^Cz^)  to  the  function 
Fij(Xj)  of  Figure  1,  for  the  choice  of  { z ° . . . z^j } shown  (here  t=6) 


Figure  2.  THE  APPROXIMATING  FUNCTION  F..(z.) 


f 


» 


When  lj  and  Uj  are  finite,  as  is  often  the  case  in  applications, 
then  it  follows  from  the  first  theorem  of  Weierstrass  that  by 
increasing  t,  and  by  judicious  choice  of  the  cut  points 

we  can  approximate  P1j(xj.)  arbitrarily  closely 
(according  to  most  of  the  standard  measures  of  "closeness"). 

In  this  way  we  have  constructed  an  approximation  to  Problem 
P-1  which  we  shall  call  P-2: 


n 

^minimize  P-,(z)  = l F^j(z^) 


where 


J=1 

z = (z1...zn) 


P-2  / 


subject  to  F^z)  = l F1  (z,  ) 


j =1  1 ^ 


Vz>  5 X W ' bi 


i . < z . < u . 
J - J - j 


i=2. 


• > Q 


i=q+l , . . . ,n 
j-1, . . . ,n. 


Program  MOGG  constructs  the  Problem  P-2  from  P-1  and  finds 
an  optimal  solution  thereof.  Interested  readers  are  referred 
to  Appendix  A,  which  discusses  this  approach  in  greater  detail 
and  which  describes  and  rigorously  justifies  the  algorithm 
employed  by  MOGG. 


C.  USER'S  GUIDE1 

This  section  provides  the  information  necessary  to  use 

Program  MOGG.  The  notation  is  from  Section  A.  We  make  the 

following  conventions.  For  any  variable  x^ , if  at  least  one 

F (x.)  is  nonlinear,  that  is,  if  it  is  not  of  the  form  a. . • 
ij  j 


JIn  this  section,  0 will  represent  "zero"  and  0 will  represent 
the  letter  "oh." 
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where  a^j  is  a constant,  then  we  will  say  that  Xj  is  a nonlinear 
variable.  Otherwise,  we  will  say  that  Xj  is  linear.  If  Xj  is 
linear,  then  Program  MOGG  assumes  £ j = 0 and  = +00.  When 
this  is  not  the  case,  then  treating  Xj  as  a nonlinear  variable 
with  1 cut  to  enforce  the  upper  bound  is  permissible.  Follow- 
ing are  the  input  specifications  for  MOGG. 

Two  types  of  input  are  required:  a user-supplied  sub- 
routine and  data  cards.  We  describe  the  subroutine  first. 


Subroutine  GETPHI 

One  component  of  input  necessary  to  use  MOGG  is 
Subroutine  GETPHI  (I,  J,  X,  F).  Called  by  MOGG,  and  given  the 
values  of  I,  J,  and  X,  GETPHI  must  set  F equal  FTT(X).  The 
value  of  X supplied  by  MOGG  will  always  equal  some  cut  zj.  The 
value  of  J will  never  correspond  to  a linear  variable.  At 
present,  no  user-supplied  read-in  capability  is  provided.  It 
is  an  elementary  matter  to  modify  MOGG  to  build  in  such  a 
capability. 


Data  Cards 


Columns 


6-10 


11-15 


16-20 


Specification  Card 
Entry 

NMROWS  - the  number  of  rows  of 
Problem  P-1,  corresponds  to  m of 
Section  A.  Note  that  this  includes 
the  objective  function. 

NUMVAR  - the  number  of  columns  of 
Problem  P-1,  corresponds  to  n. 

MAXLP  - the  maximum  number  of  calls 
to  the  linear  programming  subsection 
permitted  (100  is  a typical  choice). 

KBUB,  = 1 if  an  upper  bound  for  the 
optimal  solution  is  to  be  provided, 
otherwise  leave  blank. 


Format 
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21-25  IXPRIN,  = 1 if  the  user  wants  printed  15 

all  feasible  points  found,  otherwise 
leave  blank. 

26-30  Kl,  = 1 if  the  user  wants  all  LP  15 

solutions  printed,  otherwise 
leave  blank. 

31-35  K2,  = 1 if  the  user  wishes  to  see  the  15 

packed  LP  matrix  at  the  beginning 
of  the  run,  otherwise  leave  blank. 

36-40  K3j  = 1 will  print  LP  iteration  15 

information.  Use  for  debugging 
only — leave  blank  for  general  use. 

41-45  K4,  = 1 if  the  user  wishes  to  see  the  15 

branch  and  bound  list  after  each 
stage  is  completed. 

46-50  K5,  = 1 if  the  user  would  like  to  15 

scale  the  LP  matrix  by  dividing 
each  row  by  a power  of  2 near  the 
geometric  mean  of  the  largest  and 
smallest  (in  absolute  value)  nonzero 
entries  in  that  row. 


Upper  Bound  Card 

This  card  is  included  only  if  KBUB  = 1 on  the  Specification 
Card.  It  contains  the  user-supplied  upper  bound  (Format:  F10.6). 

Relation  Cards 

These  cards  specify  the  row  type.  Enough  cards  are 
necessary  to  allow  2*NMR0WS  columns  which  are  considered  to  be 
numbered  sequentially.  Columns  1 and  2 contain  b0  (b=blank, 
0=zero).  For  k=2,  . ..,NMROWS,  columns  £k-l  and  2k  contain 

-1  if  row  k of  the  input  problem  is  an  equality, 
bl  if  row  k is  an  inequality  (only  < is  allowed). 

Contrary  to  the  notation  used  for  Problem  P-1,  inequalities 
and  equalities  may  be  listed  in  any  order. 


These  cards  contain  the  convexity  flags.  Enough  cards  are 
necessary  to  provide  NMROWS  columns,  numbered  sequentially. 
Column  k contains  a 0 if  k=l  or  if  row  k of  the  input  problem 
represents  a nonconvex  constraint.  If  row  k represents  a 
convex  constraint,  column  k contains  a 1.  When  unsure,  the 
user  should  use  a 0. 


Bound  and  Cut  Cards 

For  each  variable,  there  is  a set  of  cards  as  follows. 
Columns  1-5  of  the  first  card  contain  the  variable  number 
(format  15)-  These  numbers  must  start  at  1 and  increase  up  to 
NUMVAR.  Columns  6-10  contain  the  value  of  the  variable  NOINC 
(format  15).  For  linear  variables,  WOINC  = 0 and  no  further 
entries  or  cards  are  required.  For  nonlinear  variables,  NOINC 
is  the  number  of  cuts  desired  for  this  variable.  NOINC  is  the 
same  as  "t"  in  Section  A.  If  NOINC  r 0,  then  columns  11-15 
must  contain  either  "AUTO."  or  "MANU."  (format  A5).  The  period 
must  appear.  If  "AUTO."  appears,  MOGG  will  automatically  make 
the  cuts.  The  next  card  must  contain  the  values  of  £j  (columns 
1-10)  and  Uj  (columns  11-20)  for  this  variable  (format  2F10.6). 
No  further  cards  are  then  needed.  If  "MANU."  appears,  then  the 
values  of  to  z ^HNC  must  appear,  in  order,  on  the  next  cards. 
Each  zj  occupies  an  F10.6  field.  As  many  cards  as  necessary 
are  to  be  used. 


Right  Hand  Side  Cards 

Enough  cards  are  required  to  provide  NMROWS  F10.6  fields. 
These  contain,  in  order,  the  right  hand  sides  of  Problem  P-1 
(the  b^).  The  first  field,  corresponding  to  the  objective 
function,  must  contain  0.0. 


Linear  Variable  Cards 


For  each  linear  variable,  in  order,  enough  cards  are 
required  to  provide  NMROWS  F10.6  fields.  The  field  contains 
the  coefficient  of  the  linear  variable  in  row  k. 

Variable  Names  Cards 

Enough  cards  are  required  to  provide  NUMVAR  A5  fields. 

These  contain,  in  order,  alphanumeric  names  for  the  variables. 

If  no  variable  names  are  desired,  then  a sufficient  number  of 
blank  cards  must  be  supplied. 


i 


Problem  Title  Card 

Finally,  one  card  must  be  provided  for  the  problem  title, 
Any  alphanumeric  expression  will  do. 


D, 


SAMPLE  PROBLEM 


’ 


s problem  i 

s discussed 

in 

Minimize 

2x? 

- 9xj  + 9x 

subject  to 

6x  j 

- l8x2 

< o 

-6x  j 

+ 18x2 

< 9 

VI 

o 

xP  x2 

< 3 

Figure  3 is  a listing  of  the  subroutine  GETPHI.  Figure  4 
reproduces  the  data  cards  for  this  problem.  Figure  5 shows  the 
MOGG  output.  This  run  took  3.7  seconds  of  CPU  time  (on  a CDC  6400 
computer ) . 


SUHKIIl'TlMr  (,F  TPHj  ( I » J,  X ,F  ) 

F = n.n 

r.Ol  ii  ( 1 IK  ,?ru,300)  , j 

In"  TFTJ.FG.i  )F».*.«X#«i-9.*x*x*9.»X 

TFIJ.Fl.  .?)F=I-2.)*X*»X»9,«X*X-9.*X 
pEtuhm 

2nr  IF  ( j.fi; , ) i F«n.*x*x 

TF  ( J.FI..71F  * (-1B.  ) *X 
PETUK" 

TF  ( J.Fr.  ) f x (.ft,  ) *x*x 
TF  ( F=  IH.*X 
pETukp 
FNIi 

Figure  3.  SAMPLE  GETPHI 
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Figure  4.  DATA  CARDS  FOR  SAMPLE  PROBLEM 


program  mogg--fiw>s  global  solutions  to  approximate  problems 
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0.  3»l'00 
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0.  0.  9,000 
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1 1 7 

2 8 1* 

STARTING  TO  iTERATt 

stage. problem  lomer  bound  upper  bouno  branching  vaRIablf 


o.n 

-3.667 

-2.667 

2 

done 

blb. 

mith  This 
-3.607 

STAGE 

, BUB*  -2.667 

, BRANCHING  on  prorlem 

0.0. 

VARIABLE 

NUMRER 

2 

1.1 

LB  GT  BUB 

1 .7 

-3.500 

-2.000 

1 

1 • 

-3.5o0 

-2.000 

1 

done 

BI.Bs 

MlTH  ThIb 
-3.60ft 

STAGE 

. RUB,  -2.667 

, BRANCHING  on  prorlem 

1.3* 

variable 

NUMBER 

1 

?•! 

-2.857 

-2.857 

0 

?.? 

LB  GT  BUR 

done 

HI.Ba 

MlTM  ThI6 
“3.60ft 

STAGE 

. nU0«  -2*857 

• BRANCHING  on  problem 

1*2* 

variable 

number 

1 

3.1 

-2.857 

-2.857 

0 

3.?  LB  GT  BUB 

--•-SAMPLE  problem 

objective  function  at  optimum  -2.B57 
variable  VALUES  at  optimum— 

X]  X? 

t.TU  1.000 

Figure  5.  MOGG  SAMPLE  OUTPUT 
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Note  some  of  the  features  of  the  output.  The  KLO,  KRO 
columns  display  the  limits  of  the  "k-sets"  of  Appendix  A,  which 
are  stored  as  a single  variable  array.  For  computational  pur- 
poses, linear  variables  are  assigned  a k-set  in  which  KLO  equals 
KRO.  MOGG  prints  "STARTING  TO  ITERATE"  after  completing  its 
data  storage  routines,  and  begins  the  branch  and  bound  proce- 
dure. Problems  are  numbered  by  their  stage  and  their  position 
in  that  stage.  After  completion  of  each  stage,  a best  lower 
bound  (BLB)  and  a best  upper  bound  (BUB),  if  any,  are  displayed. 

If  no  best  upper  bound  is  found,  BUB  will  be  set  equal  to  1.E70. 

If  no  upper  bound  is  found  for  an  individual  problem,  the  word 
NONE  will  appear.  In  Problem  1.1,  LB  GT  BUB  indicates  that  the 
lower  bound  for  that  branch  is  greater  than  the  best  upper 
bound  presently  known,  so  that  no  further  investigations  along 
that  branch  will  be  pursued.  Problem  2.1  displays  "0"  as  the 
branching  variable  to  indicate  a terminal  node  of  the  branch 
and  bound  tree. 

Additional  information  can  be  requested  on  the  specifi- 
cation card.  Most  of  the  resulting  displays  are  self  explanatory, 
however,  the  user  should  be  aware  of  the  following: 

• When  Kl=l,  the  LP  solution  will  be  printed  in 
"packed"  form  so  that  basic  variables  which 
are  equal  to  zero  will  be  omitted. 

• When  K2=l,  the  packed  (zeros  omitted)  matrix 
will  be  printed  by  columns  going  across  the 
page,  with  the  row  number  beneath  the  entry. 

An  identity  matrix  is  annexed  to  the  left 

of  the  structural  matrix. 

• When  K3=i,  the  user  should  refer  to  Appendix  B 
for  an  explanation  of  the  LP  iteration  printout. 

• When  K4=l,  the  column  beneath  "FLAG"  contains 
the  pointer  used  to  divide  the  k-sets  (x^  in 
Appendix  A). 
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E.  ON  THE  ALGORITHM 


Appendix  A contains  a thorough  description  of  the 
algorithm.  Figure  6 is  a flowchart  representing  the  MOGG 
implementation  of  this  algorithm  using  some  notation  from 
Appendix  A.  The  variable  NOLEFT  is  the  number  of  problems 
left  to  solve  in  any  given  stage. 

The  linear  programming  code  used  by  MOGG  is  described  in 
Appendix  B,  and  listed  in  Appendix  C.  It  was  chosen  for  its 
numerical  stability,  an  important  consideration  when  trying 
to  solve  "real  world"  problems. 

F.  ERROR  EXITS 

MOGG  makes  numerous  diagnostic  checks  throughout  its 
operation  and,  under  some  circumstances,  will  terminate.  When 
this  happens,  a self-explanatory  diagnostic  message  will  be 
printed  along  with  a reference  to  the  region  of  the  code  where 
the  error  occurred. 

G.  VARIABLES  AND  TOLERANCES 

These  common  blocks  provide  interroutine  communication  for 
MOGG.  Block  /FIRST/  contains  mostly  main  program  variables, 
while  /WORK/  and  /BLOCK/  are  primarily  for  the  use  of  the 
linear  programming  subsection.  Among  the  important  variables 
are  the  following  (see  Section  A and  Appendix  A for  terminology): 


KLO(I),  KRO(I): 

These  define  the  lower  and  upper 
boundaries  of  variable  I’s 
original  k-set. 

W: 

This  array  is  used  by  LINPRG  to 
return  optimal  LP  solutions  . 

LFLG : 

LINPRG  uses  this  to  indicate  infeasi- 
bility of  a subproblem. 

CUTS: 

This  array  stores  all  the  cuts  Zj  . 
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Figure  6.  MOGG  LOGIC 


Start 


Read  in  data 
Set  up  problem 
NOLEFT=l 

BUB=°°  unless  provided 


Solve  next  problem  in  this 
stage.  Set  LB  for  this  prob- 
lem equal  to  optimal  objec- 
tive function  value,  if  any 


NOLEFT  = NOLEFT  - 


Put  problem  on  list,  along 
with  LB  and  information  for 
tracing  back  up  branch  and 
bound  tree 


Create  0 


Compute  branching  variable,  if  any, 
for  this  problem  and  store  it  on  list 


Is  0 a feasible  point  'V  NO 

of  the  approximate  problem?/ 


Compute  objective  function  value 
of  F.  If  less  than  BUB,  replace 
BUB  bv  this  value 


DOES  NOLEFT  = 0? 


(continued  on 
(next  page) 


p- 


Figure  6.  MOGG  LOGIC  (con't) 


» 


» 


► 


I 


♦ 


ZLSTNO  >> 
ZLSTPA  I 
LSTKL 
LSTKR 
ZLSLB 
IBRVR 
FLAG 


A,  IA: 

B: 


Seven  arrays  that  constitute  the  list 
representing  the  branch  and  bound  tree. 
ZLSTNO  stores  stage  and  problem  numbers, 
ZLSTPA  stores  the  number  of  the  immediate 
predecessor  of  each  problem,  LSTKL  and 
LSTKR  are  the  lower  and  upper  boundaries 
of  the  k-sets  which  distinguish  this 
problem  (only  the  k-sets  relating  to 
the  predecessor's  branching  variable 
are  stored).  ZLSLB  is  the  objective 
function  value  computed  for  this  problem. 
IBRVR  is  the  branching  variable  for  this 
problem  and  FLAG  is  used  to  determine 
the  new  k-sets  when  branching  on  this  node. 

These  are  used  to  store  the  packed  LP 
array . 

This  array  stores  the  right  hand  side 
values . 
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There  are  five  tolerances  specified  by  DATA  statements 
which  are  used  by  MOGG. 

BUBTOL:  Used  to  remove  insignificant  differences  when  testing 

for  the  smallest  BUB. 

FEASTL:  Margin  within  which  9 will  be  considered  feasible. 

DIFFTO:  Removes  insignificant  differences  when  choosing  a 

branching  variable. 

DONTOL:  When  BUB  and  BLB  are  within  DONTOL  of  each  other,  the 

problem  will  be  considered  solved. 

CUTTOL:  Used  to  determine  when  the  partition  indicator  (FLAG) 

for  a k-set  falls  on  a cut. 

All  other  tolerances  are  used  by  the  linear  programming 
subsection  and  should  be  changed  with  caution.  A discussion 
of  LINPRG  tolerances  is  contained  in  Appendix  B. 

The  arrays  are  presently  dimensioned  large  enough  to  solve 
most  problems  of  interest.  If  the  user  wishes  to  redimension 
the  arrays,  he  is  referred  to  the  COMMENT  statements  at  the 
beginning  of  the  MOGG  code  (see  Appendix  C).  Note  that  the 
variables  MAXVAR , MAXCUT,  LSTMAX , MAXROW  and  MAXA  must  be 
assigned  new  values.  At  present,  MOGG  can  handle 

100  Original  variables 
1100  Total  cuts 
100  Rows 

700  Entries  in  the  branch  and  bound  list 

5000  Nonzero  elements  in  the  packed 
linear  programming  array. 

The  MOGG  routine  has  performed  well  on  a CDC  6400  with 
60-bit  words.  If  round-off  problems  appear  when  the  code  is 
implemented  on  machines  with  smaller  words,  conversion  to 
double  precision  is  recommended. 
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1.  INTRODUCTION 

An  algorithm  for  finding  global  solutions  of  nonconvex 
separable  problems  was  developed  by  Falk  and  Soland  [3]  and  Soland 
[8].  The  method  is  based  on  the  branch  and  bound  philosophy  and 
yields  a (generally  infinite)  sequence  of  points  whose  cluster 
points  are  global  solutions  of  the  problem.  The  implementation 
of  the  method  is  severely  limited  by  the  necessity  of  computing 
convex  envelopes  [4]  of  the  functions  involved  although  a number 
of  applications  of  the  method  have  been  made  (e.g.,  [5],  [9]). 

These  applications  were  possible  because  of  the  special  structure 
of  the  functions  involved  (e.g.,  concave  or  piecewise  linear). 

The  traditional  method  for  treating  separable  problems 
involves  calculating  piecewise  linear  approximations  of  the  func- 
tions defining  the  problem  and  applying  a modification  of  the 
simplex  method  to  the  resulting  problem  (see,  e.g..  Miller  [7]). 

The  modification  amounts  to  a restriction  on  the  usual  manner  of 
selecting  variables  to  exchange  roles  (basic  to  nonbasic  and  vice 
versa)  and  will  yield  a local  but  not  necessarily  a global  solu- 
tion of  the  approximating  problem. 

In  this  paper  we  present  a method  that  will  yield  a global 
solution  of  the  approximating  problem  referred  to  above.  The  method 
is  similar  to  the  Falk-Soland  algorithm  but  takes  advantage  of  the 
special  structure  of  the  resulting  approximate  problem  and  employs 
the  branch  and  bound  philosophy  to  set  up  and  monitor  the  solutions 
of  a finite  sequence  of  linear  subproblems. 


* 
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Recently  Beale  and  Tomlin  [1]  announced  that  they  have  developed 
a similar  algorithm  which  they  have  incorporated  into  their  UMPIRE 
mathematical  programming  system  [10],  oasic  idea  of  their  method 

is  the  same  as  that  of  the  algorithm  detailed  herein  although  their 
rules  for  selecting  branching  nodes  and  branching  variables  are  dif- 
ferent, being  developed  from  an  integer  programming  point  of  view 
while  ours  are  modifications  of  the  rules  developed  in  the  Falk-Soland 
method  [3]  and  its  extension  by  Soland  [8]. 

The  problem  which  we  address  has  the  form 

minimize  F (x) 
o 

subject  to  F^ (x)  4 b i = l,...,m 

i i x i L 

where  i,  and  L are  finite  lower  and  upper  bounds  respectively  on  x . 
We  assume  that  each  F^  (i=0,l,2, . . . ,m)  is  separable,  i.e., 

n 

F . (x)  = E F (x.)  i = 0,l,...,m 

1 J-l  i3  3 


and  that  each  F. . is  continuous.  As  extension  to  the  case  where 
ij 


is  piecewise  continuous  is  covered  in  Section  5. 


In  Section  2 we  define  the  approximating  problem  of  problem  Q 
and  construct  the  problem- obtained  by  replacing  each  of  the  functions 
involved  by  their  convex  envelopes.  A related  problem  is  simultaneously 
introduced  and  shown  to  give  a sharper  underestimate  of  the  optimal 
value  of  the  approximating  problem  than  does  the  convex  envelope  prob- 
lem. It  is  this  related  problem  which  the  branch  and  bound  procedure 
solves  first  to  get  estimates  on  the  optimal  value  of  the  approximating 
problem  and  to  set  up  new  problems  if  the  estimates  do  not  yield  a 
global  solution. 
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A detailed  analysis  of  the  complete  method  is  given  in  Section  3 
and  an  example  follows  in  Section  4.  Some  computational  considerations 
are  given  in  Section  5. 


2.  THE  APPROXIMATING  PROBLEM  AND  CONVEX  ENVELOPES 


The  approximating  problem  of  the  original  problem  Q is  obtained 
by  replacing  each  function  by  a piecewise  linear  approximation 

over  the  interval  . One  common  method  (see,  e.g.,  [7])  that 

is  employed  involves  selecting  p.  + 1 grid  points  y.  , . . . ,y  in 

3 J ® jPj 

Lj  ] where  y^Q  = and  y^  = L^  and  using  convex  combinations 

of  the  numbers  F^j(y^)  and  F^  (y^  as  approximations  to  the 

values  of  F (x  ) over  the  subinterval  [y  ,y  ,, , ] . Figure  1 
ij  j 3 k j > *+1 

illustrates  this  type  of  approximation. 


Figure  1.  PIECEWISE  LINEAR  APPROXIMATIONS 


Mathematically,  we  obtain  this  approximation  by  setting 
Fij(Xj)  = fij(0j)  where 


where  K.  = {0,1,, 
J 


j 

i . . ,p  . } , 6 . = (9  . 0 . ) if 

J i J°  * JPj 


ktK . 
J 


Vjk = xj 


(2.1) 


(2.2) 


i e.,  = l 

keK.  Jk 


(2.3) 


9.,  > 0 
jk  “ 


keK 


j 


(2.4) 


and  if  we  add  the  further  restriction  that  at  most  two  of  the  weights 

{0.,  : keK.}  are  nonzero,  and  if  two  are  nonzero,  then  these  must 

jk  J 

correspond  to  adjacent  grid  points.  This  last  restriction  is  necessary 

since  without  it  one  may  obtain  any  point  in  the  convex  hull  of  the  set 

{ (y . ,F. . (y,  )),... (y.  ,F..(y.  ))}  which  lies  on  the  vertical  line 

jo  ij  jo  J.Pj  1J  J,Pj 

passing  through  x^  . As  we  shall  have  occasion  to  refer  to  this 
restriction  later,  we  make  the  following  definition. 


The  Adjacent  Weights  Restriction  (AWR) : Let  KCK.  be  a set  of 

consecutive  integers.  The  set  of  numbers  {0.,  : keK}  satisfies  the 

Jk 

adjacent  weights  restriction  if  at  most  two  of  these  numbers  are  nonzero, 

and  if  0 , 0 > 0 then  either  s=t-l  or  s = t + 1 . 

js  j t 

We  shall  use  the  symbol  f..(0.)  to  denote  the  function  defined 

ij  J 

by  expression  (2.1)  and  the  symbol  f„(x.)  to  denote  the  function 

f^j(0j)  constrained  by  relations  (2.2),  (2.3),  (2.4)  and  the  AWR. 

Thus  denotes  a linear  function  of  the  variables 

0.  ,0. 0.  while  f..(x.)  denotes  a piecewise  linear  function 
J°  J1  JPj  ^-J  J 
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of  the  single  variable  such  as  that  illustrated  in  Figure  1. 


Likewise 


f (0)  = £ f (0  ) 

j = l 3 3 


f (x)  = £ f (x  ) 

J-l  3 3 

for  i = 0,1,2, ... ,m  . 

By  replacing  each  ^ij^xj^  by  its  piecewise  linear  approximation 
f„(x^)  , we  obtain  the  following  approximate  problem 


minimize  f (0)  = £ £ Q.,F.(y  ) 

0 j-l  keK.  °3 

J 


problem  P 


subject  to  f.(0)  * £ £ 0.,F..(y.,)  < b, 

1 j=l  kcK.  3k  13  3 

J 

£ 0.,  = 1 

keK.  3k 


(i  = 1, . . . ,m) 


(j  ~ 1 f • • • f n) 


Vi0 


(j  = !» # • • |H  ! keK^ ) 


I {0.,  : KeK.}  satisfies  AWR  (j  = l#...,n) 

[ jk  j 

Here  0 = (O^O^  . . . ;0n>  = (01O»  • • • •91pi;e20’  * 1 * ’ * * * ;0n0 0npn*  * 

The  solution  value  of  this  problem  is  offered  as  an  approximation  to  the 

solution  value  of  the  original  problem,  problem  Q.  The  solution  point 

£ 

0 of  problem  P yields  an  approximation  to  the  solution  of  problem  Q 
via  the  relations  (2.2),  i.e., 


J 


j — 1 , • • * ,n 


Problem  P is  the  usual  problem  that  is  addressed  when  seeking 
solutions  of  separable  programs  (see,  e.g.,  [7]).  The  method  of 
"solution"  involves  generating  a basic  feasible  solution  of  the  linear 
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program  associated  with  problem  P that  satisfies  the  AWR.  A 
modification  of  the  simplex  method  is  then  used  to  sequentially 
change  the  basis  until  a local  solution  of  problem  P is  obtained. 
This  modification  amounts  to  a restricted  basis  entry  rule  which 
insures  that  the  AWR  are  always  satisfied  by  the  basic  feasible 
solution  associated  with  each  stage  of  the  simplex  method.  Thus  the 
only  nonbasic  variables  9^  that  may  enter  the  basis  at  a given 

iteration  are  neighbors  of  existing  basic  variables.  If  such  a 
variable  is  chosen  to  enter  the  basis,  the  outgoing  basic  variable 
must  be  chosen  so  that  the  new  basic  feasible  solution  satisfies 
the  AWR.  It  may  be  shown  that  this  method  will  yield  a local 
solution  of  problem  P,  so  that  if  problem  P is  convex,  the  solution 
will  be  a global  solution.  In  particular,  if  problem  Q is  convex, 
then  so  is  P and  a global  solution  is  assured. 


In  this  paper  we  are  concerned  with  a method  that  will 
produce  global  solutions  of  problem  P.  The  method  may  be  considered 
a specialization  of  the  method  of  Falk  and  Soland  [3]  and  the  exten- 
sion described  by  Soland  [8],  In  this  method  it  is  necessary  to 
compute  "convex  envelopes"  of  all  functions  involved  in  the  problem 
description  over  appropriate  intervals.  A number  of  convex  sub- 
problems are  then  set  up  and  solved  with  the  branch  and  bound  philos- 
ophy monitoring  the  solution  values  of  these  problems  and  guiding 
the  creation  of  new  subproblems.  The  convex  envelope  of  a function 
of  a single  variable  over  an  interval  L ^ ] is  that  convex 

function  f . . defined  over  [i.,  L.l  such  that,  if  d. . is  any 
ij  1 j ’13  7 


convex  function  on  [JL,  ] 
point  in  L ] , then  d„ 


which  underestimates 
also  underestimates 


at  every 
over 


[ly  ] . Roughly,  the  convex  envelope  of  a function  is  the  highest 

convex  function  which  underestimates  that  function  over  the  appropriate 
interval.  Alternate  and  more  general  definitions  and  relations  con- 
cerning convex  envelopes  are  found  in  [4], 
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We  are  interested  in  determining  the  convex  envelope  of  the 
piecewise  linear  functions  f^(Xj)  defined  by  the  relations  (2.1) 

through  (2.4)  together  with  the  AWR.  It  is  clear  geometrically,  and 
not  difficult  to  show  analytically,  that  the  convex  envelope  of  this 
function  over  [i  , ] is  the  function  f ) : 


f,^(x  ) = min 
11  1 

keK^1^  (yjk) 

(2. 

s.  t. 

keKj  = "i 

(2. 

keK>  = 1 

1 

(2. 

0,.  > 0,  keK . . 
jk  - j 

(2. 

Note  that  we  do  not  impose  the  AWR  on  the  definition  of  f..(x.) 

il  1 

• 

We  illustrate  this  definition  in 

Figure  2 which  may  be  compared 

to 

Figure  1. 


Figure  2 . 


CONVEX  ENVELOPES 


Thus  the  calculation  of  f^Cx^)  at  a given  point  involves 

the  solution  of  a linear  program.  The  first  subproblem  addressed  by 
the  method  described  in  [8]  would  be 

n 

min  f£(x)  “ 2 fQ^(Xj) 


subject  to  f (x)  = I f.  . (x  ) <_  b 

j=l  3 J 


(i  = 1, . . . ,m) 


1 < x < L . 


This  is  a convex  program  whose  solution  value  serves  as  an  underestimate 
of  the  solution  value  of  problem  P.  Because  of  the  piecewise  linear 

c 

nature  of  the  functions  f^  , it  is  possible  to  convert  this  problem 

to  a linear  program.  This  approach,  however,  involves  explicitly 

c 

calculating  the  functions  f for  each  i and  j . Moreover,  it 

would  be  necessary  to  do  this  for  a number  of  problems  of  the  above 
form.  We  may  avoid  these  calculations  by  considering  the  related 
linear  program: 


min  f (0)  = 
« o 


1 Jsubject  to  f , (0) 


Z E 9ikFoi(yik) 
j=l  keK.  3 1 

3 


j=l  keK. 


VW  =bi  (i  = 1 m) 


(j  = !»•••> tt) 


V-° 


(j  = 1 n ; keK  J 


Note  that  problem  P is  similar  to  problem  P except  that  the  AWR  are  not 
present.  Moreover,  given  a feasible  point  0°  of  problem  P,  it  follows 
that  the  point  x°  defined  by 


6W 


is  feasible  for  the  convex  envelope  problem  by  virtue  of  the 
inequality 

fUC(«?  i Ve°>  • 

It  is,  however,  possible  that  the  convex  envelope  problem  has  feasible 
points  x for  which  there  is  no  feasible  6 satisfying  the  above 
expression.  For  if  x is  feasible  to  the  convex  envelope  problem, 
for  each  i = l,...,m  there  must  be  a vector  i which  satisfies 
conditions  (2.6),  (2.7),  and  (2.8)  together  with  the  conditions 
f^(6)  4 b^  . This,  in  itself,  does  not  imply  the  existence  of  a 
single  vector  which  satisfies  all  of  these  conditions.  On  the  other 
hand,  any  point  feasible  to  problem  P is  also  feasible  to  so  that 
the  solution  value  of  P^  offers  a valid  lower  bound  on  the  solution 
value  of  P. 

3.  THE  BRANCH  AND  BOUND  ALGORITHM 

In  this  section  we  present  an  algorithm  to  calculate  the  global 
solution  of  problem  P which  is  based  on  the  branch  and  bound  philosophy 
(see,  e.g.,  [6]).  The  algorithm  considers  subsets  of  a linear  polyhedron 
containing  the  feasible  region  F(P)  of  problem  P.  A lower  bound  on 
the  optimal  value  of  problem  P is  found  by  minimizing  fo(0)  over  each 

of  these  subsets  and  selecting  the  smallest  of  these.  A check  for 
solution  is  made  which,  if  successful,  yields  a global  solution  of  P. 

If  the  check  fails,  the  subset  corresponding  to  the  smallest  lower  bound 
is  further  subdivided  into  either  two  or  three  new  linear  polyhedra  and 
the  process  continues  as  before  with  new  and  sharper  bounds  being  deter- 
mined. The  process  is  finite  and  terminates  with  a global  solution  of  P. 

As  is  customary  with  branch  and  bound  procedures,  the  algorithm 

is  described  in  terms  of  a branch  and  bound  tree.  (See  Figure  6 for  an 

example.)  The  nodes  of  the  tree  will  be  identified  with  the  symbols 
12  3 i 

N , N , N ,...  and  each  node  N will  correspond  to  a linear  subproblem 
P*  of  problem  P.  It  is  convenient  to  also  use  the  notion  of  a "stage." 
The  first  stage  of  the  method  consists  of  problem  P^  (or  node  N^)  and 
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its  solution.  The  second  stage  of  the  algorithm  consists  of  problems 

P*  together  with  either  2 or  3 new  subproblems  created  from  problem  P*. 
A new  stage  is  created  when  a previously  solved  subproblem  is  chosen 
for  branching  and  new  subproblems  are  formed.  For  example,  the  tree  of 

Figure  6 illustrates  that  8 subproblems  were  formed  in  4 stages.  The 

1 12  3 

first  stage  contains  node  N ; the  second  contains  nodes  N , N t M 

4 

and  N ; the  third  stage  contains  these  nodes  and  the  new  nodes 

5 6 18 

N and  N , and  the  fourth  stage  contains  nodes  N through  N 

With  each  node  NC  there  is  associated  a lineai  program  of  the 


minimize  f (0)  * E £ 0..  F . (y..) 

° j-1  keKj  ik  0j 


problem  P 


subject  to  fi(6)  = £ E eiicFii^yik^  - bi  (i*l#...,m) 

j=l  keK  J 3 J 

J i 


1 61k  = 1 

keK.  J 


(j-l,...,n) 


v° 


(j=it • • • ksKj ) 

(j=l,..., n;  k^kj) 


where  the  sets  0=1,..., n)  are  subsets  of  consecutive  integers 

of  the  sets  K^  . Note  that  each  problem  is  a linear  program  and  that 

these  problems  differ  only  in  the  constraints  0.,  = 0 0=1,..., n;  MK*)  . 

J R J 


Problem  P^  has  K^  = K^  (j=l,...,n)  so  that  problem  resembles  problem  P 

except  that  problem  P^  does  not  have  the  AWR  imposed  on  it.  Let  F(PC) 

denote  the  feasible  region  of  problem  PC  and  F(P)  denote  the  feasible 
region  of  problem  P.  Note  that 


F(P)C!F(P  ) 


(3.1) 
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p ** 


and  that  F(P  ) is  a linear  polydedron  whereas,  in  general,  F(P)  is 
not  even  a convex  set.  Assuming  F(P)  + 0 , problem  P"*-  will  have 
at  least  one  minimizing  point  01  . In  general,  let  0C  denote  a 
solution  of  problem  P1",  if  one  exists,  and  set 

>t 


LB(t)  = { 


fo(0t) 


+» 


if  9 exists 
otherwise . 


It  follows  that 

LB(t)  < min  {f  (9)  : 0eF(Pt)  H F (P) } . (3.2) 

— o 

* 

It  is  sometimes  possible  to  obtain  an  upper  bound  on  f (0  ) 

t 'V 

from  problem  P . In  fact,  if  9 is  any  feasible 

'V, 

point  to  problem  P,  the  number  f q C 0 ) will  be  an  upper  bound  on 
* t 

f o ( 0 ) . Using  the  vector  9 (assuming  it  exists)  we  may,  at  little 

computational  expense,  attempt  to  construct  a vector  9 which  is 
feasible  to  problem  P according  to  the  following  rule: 

Compute  the  vector  x*"  using  the  relationship 

*3  ' ^ "jVjk  <J'1 n)  • 

j 

We  then  compute  a vector  0 which  satisfies  the  AWR  and  the  relationship 
^ ■ kcK  <J-1 n>  ■ 


This  computation  is  straightforward  since  each  x^  must  be  in  some 

interval  [y  , ,,  y.  , , . ] and  hence  may  be  expressed  as  a convex 
j , fc  j , k 

combination  of  the  two  adjacent  points  y , , and  y.  , . If  this 

j , k j , k +1 

vector  also  satisfies  the  constraints  f^(9)  £ b^  (i=l m)  , 

the  number  fQ(0  ) serves  as  an  upper  bound  on  fQ(0  ) • We  define 

A-ll 


the  quantity 


UB(t) 


&<**> 


if  is  feasible  to  P 

otherwise 


so  that 

f (0*)  < UB(t) 
o 


serves  as  a complementary  inequality  to  (3. 2). 


(3.3) 


In  general,  the  A-th  stage  of  the  algorithm  consists  of  problems 

PL  together  with  their  solutions  0^,...,0L  (if  they  exist) 

and  the  quantities  LB(1),  UB(1) , . . . ,LB(L) , UB(L)  . A node  (or  equiva- 
lently, a problem)  from  which  no  branching  has  yet  taken  place  (from 
which  no  new  problems  have  been  created)  is  termed  an  intermediate  node 
(intermediate  problem) . The  set  of  all  intermediate  problems  at  stage  A 
is  denoted  by  1(A)  . At  stage  one,  1(A)  = {1}  , and,  if  three  new 
problems  are  created  to  form  stage  two,  1(2)  = {2,3,4}  . 

The  algorithm  is  to  be  constructed  in  such  a way  that 

*(P)  C U F(Pt)  • (3.4) 

tel(A) 

I ) 

We  define  the  quantities 

BLB(A)  = min  {LB(t)  } 
tel (A) 

and 

BUB (A)  = min  {UB(t)}  . 

t — 1 , . . . ,L 

Then  (3.2),  (3.3)  and  (3.4)  imply  that 

BLB(A)  < f (0*)  < BUB (A)  . (3.5) 

o =■ 

This  is  the  basic  inequality  which  signals  the  completion  of  the 
algorithm  when  equality  is  attained  throughout.  We  will  show  that  our 
method  of  branching  (creating  new  problems)  sequentially  sharpens  (3.5) 
stage  by  stage  and  will  produce  equality  in  a finite  number  of  stages. 


1 


J 
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Check  for  Solution;  If  BLB(Jl)  = BUB(Ji)  at  the  Jl-th  stage,  an 
optimal  solution  of  problem  P is  0t  where  UB(t)  = fo(0t)  = BUB(i)  . 

If  BLB(Jl)  < BUB(£)  we  must  choose  a node  NC  for  branching, 

i.e.,  a problem  PC  to  create  new  problems  which  will  sharpen  the  bounds 
in  (3.5).  We  shall  use  the  notion  that  the  numbers  LB(t)  represent 

approximations  to  the  quantities  min  (fQ(0)  : 0eF(Pt)  DF(P)  < , Since 

we  are  interested  in  determining  min  (fQ(0)  : 0eF(P)l  , we  choose 

the  smallest  of  the  numbers  LB(t)  to  determine  PC,  the  problem  most 
likely  to  generate  a global  solution  of  P. 

T 

Choice  of  Branching  Node:  Choose  an  intermediate  node  N for  further 
branching  where  LB(T)  = BLB(i)  . 

Actually  the  algorithm  will  converge  if  any  intermediate  node 
is  selected  for  further  branching  and  it  is  sometimes  convenient  from 
a computational  point  of  view  to  use  a different  rule  for  branching. 

A common  alternative  is  to  select  that  problem  which  has  been  solved 

\ 

last  for  further  analysis,  since  the  data  defining  that  problem  are 
on  hand  and  data  needed  for  the  new  problem  are  very  similar.  This 
alleviates  the  bookkeeping  involved  and  tends  to  minimize  the  number 
of  times  a particular  branch  in  the  tree  is  revisited.  On  the  other 
hand,  the  tree  tends  to  grow  larger  than  the  tree  our  rule  would  grow 
and  would  not  be  efficient  if  the  total  time  required  is  largely  a 
function  of  the  time  required  to  solve  the  subproblems.  In  our  applica- 
tion, the  amount  of  data  required  to  distinguish  one  problem  from 
another  is  minimal  so  that  this  should  not  be  a factor. 

T 

Having  selected  node  N for  branching  at  stage  l , we  create 
new  subproblems  by  choosing  a branching  variable  0j  (or,  equivalently, 

T 

xT)  and  partitioning  the  set  K into  subsets  of  consecutive  integers. 

J 

The  rule  for  selecting  follows. 
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Choice  of  a Branching  Variable:  Compute  each  of  the  differences 


r\j  T T 

‘ <vwv 


(3.6) 


for  i ■ 0,1,..., m and  j = l,...,n  . Select  J which  corresponds 
to  the  largest  of  these  differences. 

If  all  of  these  differences  were  nonpositive,  upon  s imnur.g 
over  j for  each  i = 0,l m we  obtain 


Thus 


and 


f.(eT)  - f.(eT)  <_  0 


f (eT)  < f (eT)  = blbU) 

O =“0 


f±(eT)  4 f.(eT)  < i,± 


(i=  1, . . . ,m) 


VT 


Since  0 satisfies  the  AtfR*  we  see  chat  0 eF(P)  so  that 

,*TS 


BUB  U)  < f (9  ) < BLBU) 

= o *“ 


that  is,  0 must  have  been  a global  solution  of  problem  P,  contradicting 
our  previous  assumption.  Thus,  unless  we  are  at  a solution,  at  least 
one  of  the  differences  (3.6)  is  positive  and  we  choose  J corresponding 
to  the  largest  of  these  quantities. 

This  rule  for  selecting  a branching  variable  is  analogous  to  the 
rule  suggested  in  [3]  and  [8],  Since,  at  a solution,  all  differences 
(3.6)  will  be  nonpositive,  we  are  selecting  a variable  corresponding  to 
the  worst  violation  of  this  criterion. 


Note  that  not  all  differences  (3.6)  need  be  calculated  at  every 
stage  since  some  will  automatically  be  zero.  If  the  set  {9 


jk 


T 

>eV 


satisfies  the  AWR,  for  some  j then  0 


j 


*5 


so  that  all  of  the 


corresponding  differences  (3„6)  for  i = 0,...,m  are  zero.  Moreover, 
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if  is  a convex  function,  the  piecewise  linear  approximation 

f u > of  it  (equations  (2.1)  through  (2.4)  and  the  AWR)  will  also 
be  convex.  If  we  denote  this  approximation  by  j (Xj ) we  have 


£ ^44(y..)  = z 


jk‘ij 


jk  ij  wjk' 


a, 

J 


(9j  satisfies  AWR) 


Fij(keK.8jkyjk) 


= Z 9^  F (y  ) 
keKj  1J  Jk 


so  that  the  corresponding  differences  (3.6)  are  automatically  nonpositive. 
Incidentally,  this  also  proves  that  the  algorithm  yields  a global  solu- 
tion of  a convex  program  in  a single  stage. 

Having  selected  variable  J for  branching,  we  now  are  in  a 

position  to  create  the  new  problems  of  the  (£+l)-st  stage.  Let 

T T 

K - (p,p+l q,...,r)  . Note  x ^ y since  in  this  case 

J J Jp 

0^=1  while  6^  0^  = 0 and  the  difference  (3.6)  would 

Jp  J,p+1  Jr 

T T 

be  zero.  Likewise  xT  ^ yT  . Note  also  that  K.  contains  at  least 

J J L 

three  indices  for  otherwise  branching  could  not  take  place  on  this 

T 

variable.  We  may  assume  that  xj£tyja»  y jb J where  yJa  is  the 

T 

nearest  left  neighboring  division  point  of  x^  and  y^  is  the 

nearest  right  division  point.  We  do  not  exclude  the  case  where 

T T 

XJ  “ yJa  = yJb  • where  x^  falls  on  a division  point. 
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Recall  that  problems  p\...,PL  have  been  set  up  and  solved 


1 

I 


at  the  end  of  stage  i,  . 

Branching  Rule  (refer  to  Figure  3) : 

Let  K"  = {k  : keKj  and  yJk  < yJa> 

K°  = {a,b } 

K+  = {k  : keK^  and  y < y ,} 

Referring  to  the  general  definitions  of  problem  Pfc  at  the  beginning 

of  this  section,  define  a new  problem  P1"  by  setting  equal  to  one 

of  the  above  sets  if  that  set  contains  at  least  two  elements.  The 
t T T 

other  index  sets  are  unchanged  (i.e.,  (ji^J)).  In  this 

T 

manner  we  may  define  at  least  two  new  problems  (since  had  at 


T T 

least  three  points  and  y t x , y ^ x ) and  possibly  three  new 

JD  J J J 


problems.  These  problems  are  numbered  PL+\  PL+^  and  pL+^ 

(if  defined).  Note  that  if  a ^ b,  the  problem  whose  index  set 

= (a,b)  must  have  only  solutions  with  0*j  satisfying  the  AWR. 
J J 


The  various  possibilities  are  illustrated  by  example  in  Figure  3. 
Only  the.  first  possibility  yields  three  new  problems. 
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Kj  - {A, 5, 6, 7, 8} 


KJ 


'J4 


yJ5  yJ6 


yJ7  yJ8 

1 — a h 1 


{4,5,6}  {6,7}  {7,8} 


V\ 

T 

XJ 


H 1 


{4}  {4,5}  {5, 6, 7, 8} 


+— 7k- 


{4,5,6}  {6}  {6,7,8} 


tv-H 


{4, 5, 6, 7}  {7,8}  {8} 


Figure  3.  BRANCHING  ON  VARIABLE  Xj 


Beale  and  Tomlin  suggest  a different  branching  rule  wherein 
two  new  subproblems  are  defined  at  each  stage.  Using  the  above' 
notation,  they  set 


K"  = {k  : keKj 

and 

yjk  »yJb} 

K*  = {k  : keK^ 

*J 

and 

yja  = yJk} 

so  that  the  feasible  regions  of 

their  problems 

Pk+1 

. r,k+2 

and  P 

overlap  somewhat  more  than  ours 

do. 

Referring 

to  Figure  3,  their 

sets  K and  K*"  would  be  {4, 

5,6, 

7}  and  {6, 

7,8} 

in  the  first 

case  while  in  the  other  three  cases,  their  sets  would  define  the 
same  subproblems  as  we  do. 


In  the  remarks  which  follow  we  shall  assume  that  these  problems 

PL+1,  PL+2  and  PL+^  have  been  defined.  The  other  case  is  similar. 
We  first  note  that 

T L+3 

f(p)  n f(p  ) cf(p)  n ( u fcp1)) 

t-L+1 

T 

since  any  point  0 which  satisfies  the  AWR  and  is  in  the  set  F(P  ) 
must  be  in  at  least  one  of  the  sets  F(P^+^),  F(PL+2)  or  F(PL+^) 
Since  F(P)C  F(P^)  it  follows  that 

F(P)C  F(P)  n ( U F (PC) ) 
tel(2) 

i.e,,  F(P)C  LJ  F(PC)  . Continuing  in  this  fashion  we  verify 
tel (2) 

inclusion  (3.4): 

F(P)  C U FCP1)  (3.4) 

tel (£) 

Moreover,  since  any  point  in  one  of  the  sets  F(PL+1),  F(PL+2)  or 

L+3  T 

F (P  ) must  lie  in  F(P  ) we  have 

U F(PC)  C U FCP*1)  . 

tel  (A)  tel(^-l) 

T 

This  inclusion  must  be  strict  since  the  point  9 cannot  lie  in  any 
of  the  sets  F(P^  ),  F(P^  2)  or  F(P^+^)  . For  suppose 

T L+l  L-fl  T 

0 eF(P  ) and  Kj  = {p,...,a}  . Then  0Jk  = 0 for  k = a+l,...,r 

and  < yt  which  contradicts  the  assumption  that  x^efy.  , v . ] . 
JJa  r JJaJb 

These  remarks  yield 

F(P)C  U F(PC)  U F(PC)  5 ...  ^(P1)  (3.7) 

tel  (4)  tel(H-l) 
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t 
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t 


i.e. , the  sets  U F(Pfc)  are  converging  monotonically  towards 
ttUl) 

the  set  F(P)  . 

When  new  problems  are  created  for  the  (#.+l)-st  stage,  new 
lower  and  upper  bounds  are  calculated.  Note  that 

min  {LB(PL+1),  LB(PL+2),  LB(PL+3)}  >,  LB(PT) 

_ L+3 

since  F(P  ) jt-  F(P  ) . Moreover,  since  the  poii  r 3 for 

t=L+l 

T T 

which  £^(0  ) “ LB(P  ) is  not  feasible  for  the  new  problems,  it  is 

lll&ely  that  the  above  inequality  is  strict.  The  above  inequality, 
together  with  the  definitions  of  BLB(t)  and  BUB(JL)  yield 

BLB(l)  <....<,  BLB(i)  <_  f0(6*)  <,  BUB(i)  <,  . . . <,  BUB(l)  (3.8) 

so  that  the  upper  and  lower  bounds  are  converging  towards  the  optim-  j 
value  of  P . It  remains  to  show  that  the  process  converges  in  a finite 
number  of  stages. 


Theorem.  After  a finite  number  of  stages,  the  algorithm  yields  a global 
solution  of  problem  P. 


Proof.  At  each  stage  of  the  algorithm  an  index  Je{l,...,N)  is 

T 

selected  and  the  set  K is  subdivided  into  either  two  or  three  new 


sets  of  consecutive  integers  according  to  the  branching  rule.  Each  of 
these  new  sets  contains  at  least  two  integers.  Since  there  are  but  a 
finite  number  of  choices  for  J and  a finite  number  of  ways  of  subdivid- 
ing the  original  sets  into  sets  containing  at  least  two  consecutive 


integers,  the  algorithm  would  (if  it  continued)  eventually  produce  prob- 
lems whose  feasible  regions  contained  only  points  which  satisfy  AWR 


(i.e.,  eventually  F(P)  = ij  F(Pt))  . Such  problems  must  be  ioter- 

telU) 

mediate  problems  since  their  regions  cannot  be  further  decomposed,  and 
LB(t)  ■ UB(t)  . Thus  equality  must  eventually  occur  in  (3.8)  and  the 
algorithm  is  finite. 
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4.  AN  EXAMPLE 


Problem  Q: 

minimize 
subject  to 


Fq(x)  = + (-2x2+9x2~9x2) 


Fx(x)  ** 

-6xJ 

+ 18xz 

<,  9 

F2(x)  - 

6xl 

- 18x2 

i 0 

0 = X1‘  x2  - 3 


i. 


The  feasible  region  of  this  problem  is  sketched  in  Figure  4. 
There  are  local  solutions  near  the  points  (0,0.5),  (1.787,1.065)  and 
(2.738,3.000)  with  values  -2.50,  -2.97  and  -1.46  respectively.  The 
subdivision  points  are  taken  at  intervals  of  1/2  starting  at  0.  These 
values  and  the  values  of  functions  F^  at  these  points  are  displayed 


in  Table  1 and  the  results  of  linear  approximations  are  sketched  in 
Figure  5. 


Figure  4.  FEASIBLE  REGION  FOR  EXAMPLE 
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Each  subproblem  has  variable  0 • (0  ,0  ) - (0  ,...,0,,; 

1 2 10  16 

02O,...,026)  . The  data  provided  by  subproblems  is  given  in  Table  2 

and  the  branch  and  bound  tree  is  illustrated  in  Figure  6.  The  global 
solution  of  the  approximate  problem  is  found  to  be  the  point 

x*  « (1.714,1.000) 

with  objective  function  value  -2.857.  This  solution  is  actually  found 
at  node  6 but  not  recognized  until  problem  8 has  been  solved. 


5.  SOME  COMPUTATIONAL  CONS IDE1 ATIONS  AND  EXTENSIONS 


In  this  section  we  point  out  some  computational  aspects  of 
the  method,  some  possible  variations,  and  an  extension  to  non- 
continuous  problems. 

We  first  note  that  each  problem  P*”  contains  m constraints 

corresponding  to  the  m constraints  of  problem  Q plus  n constraints 

of  the  form  Z 0 , =*  1 . Thus  the  Generalized  Upper  Bounding  Technique 
j K 

of  Dantzig  and  Van  Slyke  [2]  may  be  used  to  advantage  here,  and 
especially  if  n is  large  compared  to  m . This  method  allows  one  to 
maintain  a basis  of  size  m x m . 

Since  each  problem  PC  is  distinguished  by  the  sets  K?  , one 

need  carry  in  memory  only  that  information  which  identifies  these 
sets,  e.g, , the  first  and  last  indices  of  the  sets.  Beale  and  Tomlin 
fl]  refer  to  these  indices  as  "flags".  The  matrix  identifying  the 
coefficients  of  the  objective  function  and  the  first  m constraints 

of  P is  common  to  all  problems  P*".  face  the  basic  solution  of  a 
problem  being  branched  from  is  not  feasible  to  the  newly  created 
problems,  it  is  not  clear  that  the  basis  of  each  problem  P1"  should 

be  carried  in  memory  along  with  the  sets  K*T  . On  the  other  hand, 

the  basic  solution  of  a problem  being  branched  from  only  fails  to  be 
feasible  to  its  descendants  by  virtue  of  one  constraint  and  hence  may 
be  useful  in  creating  basic  feasible  solutions  to  the  newly  created 
problems . 

Once  a point  6V'"  is  found  which  is  feasible  to  problem  P , 

one  could  attempt  to  produce  a feasible  solution  which  satisfies 

the  AWR  by  the  device  outlined  in  Section  3.  The  computations  necessary 

to  produce  such  a point  are  fairly  simple.  If  such  a point  8^^  may 

*v,(q  ) 

be  produced,  one  can  immediately  compute  ) and  compare  this 
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with  the  BUB(f) , updating  this  number  if  f (0^)  < BUB(i)  . 


In  such  a way  one  may  be  able  to  tighten  the  number  BUB(l)  at 


each  simplex  iteration  solving  P and  possibly  come  across  an 


* t 

optimal  solution  6 of  P during  the  solution  of  a subproblem  P . 


Of  course,  this  solution  would  not  be  recognized  as  such  until  equality 
occurs  in  (3.8). 


Finally,  we  point  out  a simple  modifcation  of  the  method  that 


will  allow  one  to  deal  with  piecewise  continuous  functions  F^  . 


In  order  to  Insure  that  problem  Q has  a solution,  we  assume  also  that 


each  F.  is  lower  semicontinuous.  The  grid  points  {y..  : keK  ; 

^ J J J 


j«l,...,n}  are  chosen  so  that  all  points  of  discontinuity  of  the 


F^’s  are  among  them.  Let  be  a point  of  discontinuity  of 


Fjj  and  set 


11m  F 


XJ  * yjK 


IJ(XJ) 


F *»  F (v  ) 

*ij  ruvyjr 


lim  F 


XJ  + yJK 


IJ(XJ) 


The  lower  semicontinuity  of  F^  at  yJR  implies  that 


F°t  £ min  (Fjj#  F*j}  • Assume,  for  the  sake  of  discussion,  that 


strict  inequality  holds,  and  define  new  indices  K , K°  and  K+ 


corresponding  to  the  quantities  F^j,  F°^  and  F^  respectively. 


These  indices  are  to  be  ordered  as 


K - 1 < K"  < K°  < K+  < K + 1 


and  corresponding  new  variables  0JR  , 0°K  and  0JR  are  defined. 
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Problem  P is  thus  redefined  with  0~  F~  4-  0°  F°  + 0+  F* 

J *v  1J  j *\  ij  j i\  u 

replacing  • 9jk  + 0jk  + 0JK  rePlacin8  eJK  and 

{...,K-1,K  ,K°,K+,K+1, . . . } replacing  , 

With  these  modifications  carried  out  at  every  point  of 
discontinuity,  the  algorithm  may  be  applied  as  before  with  no  addi- 
tional changes.  Note  that  a global  solution  of  problem  P cannot  have 

adjacent  nonzero  pairs  (0JK,  0jK)  or  (0jK,  unless  the  value 

of  F T(yT„)  is  zero,  for  otherwise  the  value  of  f could  be 

O J J K O 

decreased  by  setting  0°K  * 1 while  still  maintaining  feasibility. 
Even  if  F0j(yj^  “ 0 and  one  of  the  above  pairs  is  nonzero,  an 

equivalent  feasible  solution  may  be  found  for  which  0°  - 1 and 

J N 

which  gives  the  same  value  to  fQ(0)  • 

In  the  case  that  F°^  - F^j  continuous  from  the 

left),  one  need  only  define  two  new  variables,  say  0°  and  0*  , 

JK  J K. 

and  modify  problem  P as  above.  The  case  where  F is  right  con- 
tinuous is  similar. 
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APPENDIX  B 


A DESCRIPTION  OF  THE  LINEAR 
PROGRAMMING  SUBROUTINE  LINPRG 


Paul  F.  McCoy 


A.  INTRODUCTION 


i 


L 


I 

' 


The  subroutine  LINPRG  solves  linear  programming  problems 
by  the  standard  product  form  version  of  the  simplex  method,  as 
described  in  [1].  LINPRG  is  a slight  modification  of  the  code 
written  by  John  Tomlin  to  run  the  experiments  presented  in  [A]. 
It  was  used  again  for  the  tests  in  [2].  An  important  feature 
of  the  code  is  that  basis  reinversion  is  accomplished  by  LU 
decomposition  using  Gaussian  elimination.  The  reinversion 
algorithm  was  developed  by  Tomlin  and  is  described  in  [6].  It 
uses  a pivot  tolerance  in  choosing  the  pivot  elements  so  as  to 
compromise  the  goals  of  minimizing  the  creation  of  non-zero 
elements  and  of  pivoting  on  large  elements  to  maintain  numerical 
stability . 

B.  INTERNAL  WORKINGS  OF  LINPRG 


Reference  [3]  provides  background  reference  for  this 

j 

section . 


NOTATIONS 

NCOL  = number  of  variables  (including  structurals,  slacks  and 
art  if ic ials ) , 

NROW  = number  of  rows  (including  the  objective  row), 

x = the  (NCOL  - NROW)  vector  of  structural  variables, 

s = the  NROW  vector  of  slack  and  artificial  variables, 

c = the  (NCOL  - NROW)  vector  of  costs  (objective  function 
coefficients ) , 

A = the  [(NROW  - 1)  x (NCOL  - NROW)]  matrix  of  structural 
coefficients, 

b = the  (NROW  - 1)  vector  of  right  hand  side  values  corres- 


BASIC  PROBLEM 


minimize  cx  such  that 


Ax  I < lb 


and  x > 0 . 


ACTUAL  PROBLEM 


maximize  such  that  s >_  0,  x > 0 and 


objective  row  -*•  1 


slacks  and 
artificials 


SET-UP  PROCEDURES 


Before  calling  LINPRG,  MOGG  packs  the  constraint  coeffi- 


cients 


into  the  one-dimensional  array  A(*).  Only  non-zero  entries  are 
stored.  The  location  of  coefficients  is  maintained  by  the  row 
index  array  IA(-)  and  the  column  pointer  array  LA(*). 


A(NELEM)  = value  of  NELEMth  nonzero  coefficient, 

IA(NELEM)  = row  of  that  coefficient, 

LA(NCOL)  = the  first  element  of  A(>)  belonging  to 
column  NCOL, 

LA(NCOL  + 1)  - 1 = the  last  element  belonging  to  column  NCOL. 
The  objective  coefficients  are  placed  in  the  first  row.  The 


right  hand  side  coefficients 


are  stored  in  unpacked  form 


in  the  array  B( • ) . The  type  of  each  row  is  stored  in  array 
ISTYPE( • ) : 


ISTYPE(ROW)  = 


0 if  ROW  = 1 (objective  row) 
-1  if  equality  (=) 

1 if  inequality  (<  or  >). 


Initially  the  starting  basis  is  composed  of  the  slack  and 
artificial  variables.  On  subsequent  calls  of  LINPRG,  the  last 
basis  of  the  previous  problem  is  used  as  the  starting  basis 
with  those  variables  excluded  from  the  basis  by  MOGG  replaced 
by  the  correpsonding  slack  or  artificial.  The  basic  variables 
are  doubly  indexed  by  the  arrays  JH ( * ) and  KINBAS(-). 


JH(ROW)  = that  basic  variable  that  pivots  on  row  ROW 

KINBAS(NCOL)  = pivot  row  of  variable  NCOL  if  it  is  a basic 
variable;  0 otherwise. 


1 . Major  Subroutines 

LINPRG  uses  12  subroutines — eight  are  major,  three  are 
bookkeeping,  and  one  prints  out  the  iteration  path.  The 
eight  major  subroutines  form  the  component  parts  of  the  sim- 
plex cycle  with  LINPRG  linking  them  together.  Each  cycle 
through  the  following  flowchart  corresponds  to  one  cycle  of 
the  simplex  method  with  a basic/nonbasic  variable  interchange. 
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INVERT  starts  with  the  list  of  basic  variables  stored 
in  the  array  JH(*)*  Using  the  corresponding  coefficients 
stored  in  array  A(*),  it  calculates  the  inverse  of  the  basis 
(denoted  by  B)  using  LU  decomposition.  The  procedure  is  des- 
cribed in  detail  in  Reference  [6]. 

In  general,  the  matrix  of  basis  coefficients,  B,  is  first 
decomposed  using  Gaussian  elimination  into  the  product  of  a 
lower  triangular  matrix,  L,  and  an  upper  triangular  mitrix,  U: 


B = LU 


and  B-1  = U-1L  1 


Once  this  is  done,  a representation  of  the  basis  inverse  is 
immediate  since  the  inverse  of  a triangular  matrix  is  a simple 
rearrangement  of  the  matrix  itself.  As  an  example 

1-1 


U 


-1 


U11  u12  u13 


u22  U23 


33 


l/u1]L  0 0 


1 0 
0 1 


1 -u12/u22 


1/u 


22 


1 0 -u-^/u^ 
0 1 -u^/u^ 


0 0 1/u 


33 


and  likewise  for  L 


-1 


The  LU  decomposition  of  the  basis  is  not  unique  and  one 
wants  to  choose  that  one  which  (1)  minimizes  the  number  of 
nonzero  entries  so  that  storage  requirements  are  reduced  and 
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the  number  of  computation  in  the  BTRAN  and  FTRAN  operations  are 
minimized;  and  (2),  involves  division  (e.g.,  1/u  , l/u??, 

l/u^)  by  numbers  as  large  as  possible  to  minimize  the  growth 
of  errors  (improve  numerical  stability).  The  search  for  such 
a decomposition  is  guided  by  the  tolerance  ZTOLPV  which  will  be 
described  in  Section  3 • 

As  shown  above,  the  representation  of  the  basis  inverse 
can  be  written  as  the  product  of  elementary  column  matrices 
(often  called  eta  vectors): 


B ^ = E,  . , . E~E-, 


The  eta  vectors  are  stored  in  the  one-dimensional  array  E(*). 
The  location  of  coefficients  is  maintained  by  the  row  index 
array  IE( • ) and  the  eta  vector  pointer  array  LE( • ) . 

E(NELEM)  = value  of  NELEMtl_1  nonzero  coefficient, 

IE(NELEM)  = row  of  that  coefficient, 

LE(ETA)  = the  first  coefficient  of  E(«)  belonging  to 
the  eta  vector  ETA. 

b . FORMC  (Form  the  Cost  Row) 

FORMC  checks  to  see  if  any  variables  are  at  a negative 
value  and,  if  so,  computes  the  Phase  I objective  and  stores  it 
in  work  region  Y(-).  Otherwise,  it  stores  the  Phase  II  objec- 
tive, Y ( 1 ) =1,  in  Y(- ) . 
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c .  BTRAN  (Backward  Transformation) 


BTRAN  computes  the  ir  vector  (multipliers)  and  stores  it 
in  Y( • ) : 


Y ( tt ) = Y (Et.  . .E2E1) 


L 


objective  passed  from  FORMC 

BTRAN  is  called  the  "backward  transformation"  since  it  processes 
the  elementary  transformation  matrices  in  the  reverse  order 
in  which  they  were  created. 


d .  PRICE  (Price  Out  the  Nonbasic  Variables) 


PRICE  computes  the  reduced  cost,  dj , for  those  columns  of 


the  coefficient  matrix  eligible  to  enter  the  basis  (nonbasic 
and  not  excluded  by  MOGG): 


dj  = Y(n)A(j) 


, th 


(where  A( j ) is  the  j column  of  the  coefficient  matrix). 
PRICE  then  selects  that  column  which  will  enter  the  basis, 
JCOLP . 


e .  FTRAN  (Forward  Transformation) 


FTRAN  updates  the  column  of  coefficients  corresponding  to 
the  incoming  column,  JCOLP,  and  places  the  result  in  Y(*): 


Y = Et. . .E2E1A(JCOLP)  . 


FTRAN(l)  is  the  normal  FTRAN  described  above.  FTRAN(2)  uses 
only  the  elementary  matrices  associated  with  the  upper  tri- 
angular factor  of  B and  is  used  only  in  the  subroutine  INVERT. 
FTRAN  is  called  the  "forward  transformation"  since  it  uses 
the  elementary  transformation  matrices  in  the  order  in  which 
they  were  created. 
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f . CHUZR  (Choose  That  Row  Whose  Basic  Variable  Leaves 
the  Basis) 

CHUZR  finds  the  pivot  row,  IROWP,  using  the  ratio  tests 
described  in  [3]. 

g . UPBETA  (Update  the  Values  of  the  Current  Basic 
Variables) 

UPBETA  updates  the  current  basic  variable  values  stored 
in  array  X( • ) so  that  they  correspond  to  the  new  basis. 

h . WRETA  (Write  Eta) 

WRETA  computes  the  new  eta  vector  (elementary  matrix 
^t+l)  anc*  adds  to  representation  of  the  oasis  inverse, 

array  E( • ) : 


1 n- 


Jt+i 


:0 

nP 

Oi'-- 


where 


ni  = 


nm  1 


1/Y( IROWP) 


for  i = IROWP 

-Y( i )/Y( IROWP ) otherwise. 

(Actually,  the  divisions  are  done  only  when  Efc+1  is  used). 

2 . Bookkeeping  Subroutines 

a . SHIFTR  (Shift  Values  in  the  Work  Regions) 

LINPRG  has  two  work  region  arrays,  Y(*)  and  YTEMP(*). 
Subroutine  SHIFTR  can  shift  around  the  values  of  any  of  the 
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following  four  arrays: 


12  3 4 

B(  • ) X(-)  Y(  • ) YTEMP( • ) 

For  example,  SHIFTR(1,3)  places  the  values  of  B(*)  into  array 
Y(-),  while  SHIFTR^^)  places  the  values  in  YTEMP(-)  into 
array  Y( • ) . 


b . UNPACK  (Unpack  a Column  of  Coefficients  from  the 
Constraint  Matrix) 


Subroutine  UNPACK( JCOLP)  unpacks  the  coefficients  of  column 
JCOLP  and  places  them  in  array  Y(*). 

c • SHFTE  (Shift  Element  of  Array  E(-)) 

SHFTE  is  a bookkeeping  subroutine  used  by  INVERT.  It 
is  used  to  manipulate  the  elementary  transformation  matrices 
associated  with  the  upper  and  lower  triangular  factors  of  B. 

C.  LINPRG  OUTPUT 

Figure  1 is  an  example  of  output  generated  by  LINPRG, 
most  of  which  was  produced  by  the  subroutine  ITEROP. 

ITCOUNT  = iteration  number  (one  cycle  of  the  simplex  method 
is  an  iteration. 

I if  current  basis  is  infeasible. 

F if  current  basis  is  feasible. 

U if  solution  is  unbounded, 

blank  if  current  basis  is  optimal. 

OBJ  VALUE  = the  current  value  of  the  objective  function  (if 

STATUS  is  I,  OBJ  VALUE  is  the  sum  of  infeasibili- 
ties) . 

VECIN  = the  nonbasic  variable  coming  into  the  basis. 

VECOUT  = the  basic  variable  leaving  the  basis. 

DJ  = the  adjusted  cost  of  the  variable  coming  into  the 
basis . 

NETA  = the  number  of  eta  vectors  which  form  the  current 
representation  of  the  basis  inverse. 
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STATUS  =\ 


. 


NELEM  = the  number  of  nonzero  elements  which  form  the 
current  representation  of  the  basis  inverse. 

TIME  = 0.00,  as  the  program  timer  currently  is  not 
connected . 

IROWP  = current  pivot  row. 

X(I)  = the  adjusted  right  hand  side  on  row  IROWP. 

Y(I)  = the  current  pivot  element. 

STABILITY  COUNT  will  be  explained  in  the  next  chapter. 

Whenever  INVERT  recalculates  the  inverse  of  the  basis,  it 
prints  those  statistics  listed  in  Figure  1 under  INVERT 
STATISTICS.  These  statistics  relate  to  the  LU  factorization 
of  the  basis  and  should  be  of  little  concern  in  running  normal 
problems . 

D.  TOLERANCES  AND  OTHER  CONTROLS 

LINPRG  uses  preset  tolerances  to  reduce  the  computer 
running  time  and  accumulated  error.  These  tolerances  may 
need  to  be  adjusted  as  the  code  is  run  on  different  problems 
or  computers.  This  section  is  an  attempt  to  explain  what  these 
tolerances  do  and  how  they  should  be  adjusted. 

Solving  large  linear  programming  problems  involves  adding, 
subtracting,  multiplying,  and  dividing  many  numbers.  On  any 
digital  computer  there  are  round-off  errors  involved  in  repre- 
senting numbers  and  in  using  them  in  operations.  For  most 
programs  the  precision  of  the  computer  is  such  that  the  accumu- 
lated error  is  negligible.  Unfortunately,  most  linear  pro- 
gramming algorithms  are  designed  such  that  operations  are 
performed  on  the  results  of  operations  and,  when  this  is  done 
often  enough,  the  accumulated  error  can  grow  to  significant 
levels  even  on  precise  machines.  LINPRG  uses  the  revised 
simplex  method.  It  carries  along  a representation  of  the  inverse 

of  the  current  basis,  B-^,  which  is  a product  of  past  computations 
and  has  with  it  an  accumulated  error. 

B = ^t^t-1 ’ * ’ ^1  * 
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(Et  is  the  most  recently  added  elementary  column  matrix.)  Each 
time  the  basis  is  changed,  a new  elementary  matrix  is  added 

to  B-1  and  with  it  possibly  some  error.  At  some  point  the 

errors  may  get  out  of  hand  and  B-"*"  will  no  longer  be  a good 
approximation  to  the  inverse  of  the  basis.  Since  the  algorithm 

is  vitally  dependent  on  B_1,  it  can  then  wander  off  and  do 
ridiculous  things. 

The  accumulated  error  is  a function  of  the  number  of 
computations  and  the  size  of  the  round-off  error  involved  in 
those  computations.  In  general,  the  tolerances  allow  the  code 
to  neglect  insignificant  numbers  and,  when  choice  is  possible, 
to  perform  those  computations  with  the  smallest  round-off 
error . 

1 . Tolerances  Used  in  LINPRG 

ZTOLZE  is  the  zero  tolerance  used  throughout  the  program. 
Its  purpose  is  to  zero  out  any  "background  noise"  and  thereby 
reduce  storage  requirements  and  the  number  of  computations. 

It  should  be  slightly  larger  than  the  precision  of  the  machine — 

for  our  machine  this  is  2~ ^ « 10~^®.  If  set  too  low,  storage 
requirements  will  be  significantly  increased.  If  set  too  high, 
"good"  numbers  will  be  thrown  away  and  accuracy  reduced. 

ZTOLCR  is  the  pivot  tolerance  used  in  CHUZR.  CHUZR 
selects  the  old  basic  variable  that  leaves  the  basis  and 
thereby  the  divisor  (called  the  pivot  element)  which  is  ad- 
joined to  the  representation  of  the  new  inverse.  That  divisor 
must  have  a magnitude  greater  than  ZTOLCR.  This  keeps  the 
algorithm  from  dividing  by  small  numbers  and  thus  creating  large 
ones  which  would  increase  the  chance  of  round-off  error  in 
subsequent  computations.  If  ZTOLCR  is  set  too  high,  the 
algorithm  may  go  to  an  infeasible  basis  from  a feasible  one; 
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it  may  even  terminate  with  an  unbounded  solution  when  this 
should  not  be  the  case.  If  ZTOLCR  is  set  too  low,  errors  will 
grow  rapidly  when  the  algorithm  is  run  on  problems  which  are 
inherently  unstable. 

ZTOLPV  is  the  absolute  pivot  tolerance  used  in  the  re- 
inversion subroutine  INVERT.  It  functions  in  essentially  the 
same  way  as  ZTOLCR.  Increasing  ZTOLPV  increases  the  minimum 
size  of  the  pivot  elements  in  the  new  representation  of  the 
inverse  and  thereby  increases  the  stability.  Decreasing  ZTOLPV 
will  allow  the  representation  to  have  fewer  nonzero  elements 
and  will  decrease  the  number  of  computations  required  by  the 
algorithm.  The  tests  of  Reference  [6]  suggest  the  following 
value : 

ZTOLPV  = (102  • max|alj.  |)_1 

where  a„  are  the  coefficients  of  the  current  basis. 

ZTOLPV  and  ZTOLCR  are  related  in  that  the  revised  simplex 
part  of  the  code  hands  over  a basis  to  INVERT  which  is  non- 
singular with  respect  to  the  pivot  tolerance  ZTOLCR.  If  ZTOLPV 
is  greater  than  ZTOLCR,  then  the  reinversion  subroutine  INVERT 
may  find  that,  from  its  viewpoint,  the  basis  is  singular  and 
can  not  be  inverted.  Setting  ZTOLPV  less  than  or  equal  to 
ZTOLCR  will  avoid  this  problem. 

ZTCOST  regulates  the  tightness  of  the  terminating  test. 

If  the  minimum  adjusted  cost  is  within  ZTCOST  of  zero,  then 
the  algorithm  terminates.  It  should  be  noted  that  this  toler- 
ance does  not  affect  stability.  If  it  is  too  large,  the  solu- 
tion returned  upon  termination  may  not  be  optimal.  If  it  is 
too  small,  the  computer  time  will  become  excessive  as  back- 
ground noise  dominates. 
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No  matter  how  well  the  tolerances  are  set,  at  some  point 
the  accumulated  error  will  grow  to  significant  levels.  When 
this  happens,  the  representation  of  the  basis  inverse  should 
be  recalculated.  This  is  done  by  the  subroutine  INVERT. 
(Reinverting  is  expensive  in  terms  of  time  and  should  be  done 
only  when  necessary.  To  identify  when  it  becomes  necessary  is, 
in  itself,  a problem. ) The  accumulated  error  could  be  calcu- 
lated directly  by  computing  ||B-1B  - I||ot  . Unfortunately, 
this  would  take  about  as  much  time  as  reinverting  the  basis 
itself,  as  the  FTRAN  subroutine  would  have  to  be  called  for 

each  column  in  the  basis  in  calculating  B-1B. 

An  indirect  measure  of  the  accumulateu  error  which  takes 
relatively  little  time  to  compute  is  the  STABILITY  COUNT,  which 
appears  in  MOGG  output  for  every  call  of  LINPRG.  It  is  com- 
puted as  follows.  In  the  subroutine  PRICE  the  adjusted  cost, 
d.,  is  calculated  for  each  nonbasic  variable,  and  the  most 

t ) 

negative  is  then  chosen  to  enter  the  basis. 

adjusted  cost  = d.^BTRAN)  = c.  - c,B-'*'A(j)  . 

J J ^ 

c.  is  the  cost  for  variable  j;  c,  is  the  vector  of  the  basic 
J b 

costs;  and  A ( j ) is  column  j of  the  coefficient  matrix.  This 

is  done  by  using  the  subroutine  BTRAN  to  compute  the  multipliers. 

t r = c.  B 1 = c,  E,  . . . En  . 
b b t 1 

This  is  done  once  and  ir  is  then  applied  iteratively  to  each 
A ( j ) to  compute  the  adjusted  cost  values  in  PRICE.  tt  is 

calculated  by  multiplying  c^  and  E^  and  then  the  result  by 
Efc  ^ and  so  on.  Notice  that  the  elementary  matrices  are 
multiplied  from  left  to  right,  thus  BTRAN  is  called  the 
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"backward  transformation."  Once  PRICE  has  selected  the  non- 
basic  variable  to  enter  the  basis,  CHUZR  selects  the  variable 
to  leave  the  basis.  It  needs  the  adjusted  coefficients  of  the 
incoming  variable 

*(j)  = B-1 A( j ) = Et  ...E1A(j)  . 

This  computation  is  done  by  PTRAN  which  multiplies  the 
elementary  matrices  from  right  to  left;  thus  FT'RAN  is  called 
the  "forward  transformation."  Notice  that  with  just  one 
vector  multiplication,  dj  can  be  recalculated 

dj  (PTRAN)  = Cj  - cb'5(j)  . 

Theoretically,  matrix  multiplication  is  associative  and  thus 

d. (BTRAN)  should  equal  d.(FTRAN).  The  only  way  they  can  be 
J J 

unequal  is  if  the  "stability"  of  the  representation  of  the 
basis  inverse  has  degraded  to  the  point  where  the  accumulated 
round-off  error  generated  in  using  it  becomes  significant. 

The  STABILITY  COUNT  is  the  number  of  times 

| d . ( BTRAN ) - d. (PTRAN) | > ZTOLZE 
J J 

occurs.  Notice  the  rather  subtle  point  that  the  STABILITY 
COUNT  says  nothing  directly  about  whether  the  current  repre- 
sentation of  the  basis  is  accurate.  What  it  does  say  is  that 
if  you  multiply  a vector  by  that  representation  the  result  will 
be  affected  significantly  by  round-off  errors.  Since  the 
elementary  matrix  (which  gets  added  to  the  representation  of 
the  inverse  at  each  iteration)  is  a product  of  such  a calcu- 
lation, it  is  likely  that  it  will  also  be  in  error. 

It  has  been  found  experimentally  that  PTRAN  accumulates 
significantly  less  round-off  error  than  BTRAN.  For  this  reason 
when  mismatches  occur,  d^ (PTRAN)  is  probably  mere  accurate  than 
dj (BTRAN).  If  both  d^ (BTRAN)  and  dj(FTRAN)  are  negative,  then 
the  entering  nonbasic  variable  should  decrease  the  objective. 
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If  d.(FTRAN)  turns  out  to  be  positive,  then  we  are  probably 
J 

going  in  the  wrong  direction.  At  this  point,  LINPRG  returns 

to  PRICE  to  try  again.  If  dj(FTRAN)  turns  out  to  be  positive 
once  again,  then  the  basis  inverse  is  recalculated  by  calling 
INVERT. 

There  are  two  other  reasons  for  reinverting  the  basis. 

At  each  iteration  of  the  simplex  method,  an  elementary  matrix 
gets  added  to  the  representation  of  the  inverse.  At  some 
point  the  storage  space  will  be  exceeded  and  one  must  recalcu- 
late the  basis  inverse  representation.  The  storage  space  is 
especially  critical  for  all-in-core  codes  like  LINPRG.  The 
other  reason  for  reinverting  is  to  improve  the  running  time. 

As  the  representation  of  the  inverse  gets  larger,  it  takes 
more  computer  time  to  use  it.  At  some  point  it  will  become 
advantageous  to  expend  time  reinverting  the  basis  to  reduce 
the  size  of  the  representation  of  the  inverse  and,  thus,  the 
time  required  to  use  it  in  the  simplex  method. 

Figure  2 roughly  illustrates  how  the  computer  time  that  it 
takes  to  complete  an  iteration  of  the  simplex  method  will 
increase  as  the  size  of  the  representation  of  the  inverse  in- 
creases. Figure  3 illustrates  that  the  time  it  takes  to  rein- 
vert the  basis  will  generally  remain  constant  once  the  algorithm 
reaches  Phase  II.  The  key  question,  of  course,  is  when  should 
the  basis  be  reinverted  to  improve  the  running  time.  The  best 
solution  is  to  access  the  program  timer  and  keep  track  of  the 
time  it  takes  for  each  iteration  and  then  reinvert  according 
to  some  rule,  such  as: 

Reinvert  at  iteration  I if 
1/2 ( ITIM ( I ) - ITIM(O))  x I > INVTIM  . 

LINPRG  does  not  use  such  a rule,  since  program  timers  are 
machine  dependent  and  make  it  difficult  to  transfer  the  code 
from  machine  to  machine.  Currently,  LINPRG  reinverts  the 


I 


basis  at  least  every  50  iterations.  This  appears  to  be  a 
reasonable  approximation  to  a more  sophisticated  rule  such 
as  described  above. 

3 . Tolerance  Values 

Indications  of  tolerance  problems  are: 

1.  STABILITY  COUNT  greater  than  zero, 

2.  A large  change  in  the  objective  value  after  reinversion, 

3.  Algorithm  goes  from  a feasible  to  an  infeasible  basis, 

4.  Algorithm  takes  more  time  than  it  should  or  goes 
unbounded  when  it  should  not . 

The  principal  point  to  remember  is  that  if  the  algorithm 
does  not  behave  as  it  should,  the  tolerances  should  be  adjusted. 
If  one  wants  increased  confidence  in  the  solution,  check  to  see 
that  the  stability  count  is  low  or  at  least  that  mismatches 
do  not  occur  near  termination.  If  mismatches  occur,  adjust 
the  tolerances  and  run  again.  If  degeneracy  occurs,  perturbing 
the  right  hand  side  by  a small  amount  should  help  (this  is  called 
"epsilon  perturbation"). 

The  following  is  a list  of  values  for  the  tolerances.  The 
Orchard-Hays  column  are  values  suggested  by  Reference  [31-  The 
LINPRG  column  are  those  values  used  in  LINPRG.  The  range  val- 
ues are  the  author's  estimates  of  reasonable  upper  and  lower 
values  for  the  tolerances.  It  should  be  emphasized  that  the 
arguments  for  setting  tolerance  values  are  generally  heuristic 
and  best  values  will  vary  from  problem  to  problem. 
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Orchard-Hays 

LINPRG 

Range 

ZTOLZE 

10-12 

IQ"10 

10-18 

£ 10'8 

ZTCOST 

— 

10-10 

0.0 

< 10-2 

ZTOLPV 

10-12 

10-6 

10-10 

< 10"3 

ZTOLCR 

10"5 

lO-4 

10-10 

< 10'3 

INVFRQ 



50 

10 

< 100 

fi 
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pROGRA*  M06G(  InRLT  ,0UlRUl ,TaPE6»0UTpoT) 

•••••MaXVAr.THE  MAXIMUM  number  OF  ORtGINAL  VARIABLES*  DIMENSIONS 
KLO,KRO*K|.*K«,XCUUEU,XBESTt  VaRnmM 
•••••maxcut*  The  maximum  number  Of  internal  variables  dimensions 
w.CUTs 

.••••lstmax.  the  maximum  list  Size  dimensions  zlstno, j^st^a 
LSTKL.  LSTkR.ZLSLB  *IBRVR,FLAG 

••••*maXRO*«MAx!MUM  NUMBtR  OF  R<j"S  (rNC.  OBurCTlVE  FUNCTION), 

dimensions  icmk 
*****maXa  DIMENSIONS  a»ia 

cOmmON/FIHST/KlU( 100) *KRO ( 10U ) *<L ( 1 o ' ) *KRTloo) *XC00ED(100) *x«est 
1(100)*'*I1100)*COTS(1100)»ZlSTnO(T00)  .ZLSTRA  1 700  > * LSTKL  <r00)  , 

7 LSTKR(TOO) , ZLSLB ( 700 > »I«RVR (70o > *FL iG(7q6> IKBLI3) ,KB«(3) * 

3VARNAM (100) • PROBN  A ( 8 ) * MAX VAR  * MAxCuT * i S TMAX *«aXRO« t mAXA * 

4 NMROWC,nUmVAH,IcHK<100) »VAL*LFlG 

tUmmOn/BORkI/  B<35o>*Xl3i>O)*Y(35O),YTFMP<35o)»A<5Ooo>*E<50OO)* 

1 IA(5000),1E(5000),LA(1J02) ,LE(200?) , ICn*M u 302,2) ,K1nBA| ( i 302, , 

? JHI350) * 1ST YPE ( 35o ) , NAME ( gu ) *nTEMP<20> *CMiN,CONO,ERMAX* IFFEZ* 

3 InVfPQ.IOBJ* l^c;"P,llLH*lTCHAtiTcNT.lTRFRO.TVIN*ivOUT,jcOLP,KTNP* 

4 XSTAT,NROw,NcOl,NElEM,NETA,NLELEM,MLETA»Nv»FLEMtNGETA,NUELEM, 

F NlJFTfl  .SUmINF  ,N3 

(-Omm*N/BLOck/  ZlrLZE»^iOLPV,ZTCoST,NQMAX,NT«AX*NEMAX,URO*WMA,Q8A* 
l QFI,l.EO,QBL.ORL*QMl,QA,C)B,OC'*OE,OF.OG,QM,wTtUL,UM.WNtOO,QR,Ou,OZ 
pATA  T08J*INVfRC*ITRFMQ/1 *5o»loOO/ 
pATA  ztolZe*2Tolpv,ztcosi/i.e-Io* I.F-6»l*E-*«/ 
pATA  nmmAX,nTmaA,nemAX/35q,2uOO,5000/ 

pATa  QOO,OMA,UfciA,OFI,UtU,QBL.QPL*QM!/*MROR  , ahMATH, 4HBASI , 

1 4HFlWS,4HE0F  ,4M  *4*7  ♦ *4H  “ / 

nAT  A Q»,OB,OC*Qt,OF.OG*OM*QI ,QL,UM,0n,00,UR*oU*UZ/4HA  **HB  * 

1 4HC  , AhE  *4HF  *4HG  ,4HH  ,4Hl  »*hL  ,4HM  »4HM  , 

? 4hO  » 4hR  *4HU  *4M  Z / 

pATA  BUBTOL*FEaSTL.DlFETO*OONTOL/l.F-10* 1 »E-A, 1 .E"6* l.E-6/ 
pATA  CUTTOl/1 .E-P/ 
mAXVA«,100 
MAXCUT-1 100 
l STmaX=7oO 

mAXPOWs! 00 
wAX A*S000 
wMM«5hyANU, 
rU0»7 ,p70 

c **«*«This  section  reads  in  data 
print  t 

5 fORmaT  (6iHiPR°GRaM  mOUU„.F I NQS  gEObai  SOLUTIONS  to  APPROXIMATE  pP0 
1 pLEMS) 

o£ADin*NMR0RS*NUFV AR,MAXlP,kBUB,IXPBtN*k1*K*»K3*K*,k5 

10  FORM  a T ( 1 0 1 5 ) 

tF(krum.ne.o) Read  12, bub 
) 2 FORMAT (F10. 6) 

pHINT  15 

is  format ( i ho , i 9mphcblem  information,/) 

PRINT  20,NmRO«*s 

*0  fORm;T(1h  ,20X,110,4HKU»S) 

PRINT  ?5 , numv ar 

?5  FORMATdH  ,20X,Ii0*<)HVARlAbLtS) 

PRINT  30.MAXLP 

30  fORmaT(1H  ,20X, I ] 0,5X»26MLP  PROrLEMS  WILL  Bt  SOLVED) 

rF (kBUh.NE.o) PRINT  35»BUb 
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35 

rORMATdH  ,20X,22MUSE«"SVPPLIED  BUB 

TF<I*PfiIN.NE»0>KlNT40 

TS”,Fiu.6) 

40 

form*t<ih  ,2ox,»9hThe  ustR  requests 

1 BE  PRINTED) 

T F ( k l .nE.O)  PRlN  U5 

THAI  ALU  FEASIBLE  POINTS  FOUNO 

45 

fORmaT(1h  .20X,5cHThe  UStR  RECASTS 

1 FDj 

IF (K2.NE.0) PRINT  50 

THAI  ALU  LP  SOLUTIONS  BE  PRINT 

eO 

FORMAT (1H  ,20A,44HTHE  UStR  REQUESTS 
lF(K-».ME.O)PRlNl  5i 

THAT  THt  MATRIX  BE  PRINTED) 

51 

fORmaT(1H  ,20X,43MTHE  USER  REQUesTS 

iF<k4.nE.O)PRInT  s? 

■ P InFOHmaTION  bE  PRINTED) 

52 

fORmaT (lH  ,20Xt66MTHE  USER  R£QUESTS 
in  after  each  stage) 
iF(K5.nE.0)prInI  53 

THAT  THt  ENTIRE  LIST  BE  RRINTE 

53 

fORMaTIIH  ,20X*43hTHE  USfcR  REQUESTS 
BEAD  55, (IsTYPE i I) »I«1 *NMRQWS) 

THAT  THt  MATRIX  BE  SCALED) 

55 

fORmAT (4012) 
pRlNTG'1 

*0 

fORm;T(1H0,10HRU4  TYPE—  \ 

PRINT  65»  (lSTyP£(I) ,i«1,nmR0wS) 

65 

fORmATIIH  • 4012) 

p£An  7",  (ICHK(I) ,I»1,NMH0WS) 

70 

fORmA  T ( 80 1 1 > 

PRINT  75 

75 

FORMAT  <1Ho.17HcO('VEXI  17  FLAGS—,/) 
PRINT  «0»<lCHK<i)»I»l»NMKO*S> 

80 

FORMAT ( 1 H ,8011) 

*#***NCW  Set  up  CUTS  VtcTOR*  KLO»  ANn  KRO— 

PRINT  <50 

RO 

FORMAT (lH0,27HVARlABLt  CARDS  REpRODucEU— */ ) 

DO  100  I«1»NUMVAR 

RfAD  1 05,nOVAR,nU1nC,W0RD 

105 

FORMAT  < 15* I5»A5) 
tF(novar.ne.i>call  errU) 

print  iio,novaR,noinc,word 

1 1 0 

format 1 1 h *i5»I5»as>‘ 
if (NOINC.EU* °) 115*120 

115 

IF  (I. EQ.l)  Hb.ll7 

116 

Kl_0<  I)  *1 

KBO(I)*l 

GO  TO  100 

117 

I x*kRo ( I“1 ) ♦ 1 

rF (IX.GT.MAXCUT)CALL  err  1 2 ) 
KLOJD-KRoll-l)*1 

KBO(I)*KLO«l) 

GC  TO  100 

1 2o 

IF  ( I • EO* 1 ) 122*124 

122 

KL  o ( I > *1 

GO  TO  126 

124 

I X *KRO ( I “ 1 ) ♦ 1 

fF(IX.GT.MAXCUT»CALL  ERR(2) 
KL0(I)«KR0d-l)*l 

126 

f F < (RLO ( I) *nOInC) .GT.MAxCUT) call  ERR ( ? ) 

KCO<I)»KLO<I>*NOINC 

If  <BORD.EQi*MM)GO  TO  US 

1 1 "KLO  1 1 ) 

I 2«kRQ  < I ) 
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» MEAD  l30fCUTS(U)  ,CUTS<I2) 

1 30  FCRHAT(2Ei6.6) 

print  135,CtTS(I  U *CUTSd2) 

1 35  FORMAT ( 1h  »2GlO.*> 

IF((I2-I1).E®.1)UU  TO  100 
I*»I?-I1-1 

00  U0  ja  1 % 1 A 

t 1«0  CUTS ( 11 ♦ j) aCUTS( 1 1 > ♦ J* (CUT« ( 12) “CUTS  (111) /NO INC 

6C  TO  100 

1*5  CONTINUE 

C ***.*HERE  IF  «E  A«E  TO  Rt  A0  IN  CUTS  -ANUaLLI 

I w*KL0 ( I > 

IZ»kRO< I) 

HEAD  l50» <CoTS( J> t J*I», IZ) 

: 150  FORMAT (8F10, 6) 

PhInT  155, (CUTS (J» *J*IMt IZ) 

1 55  FORMATIlH  ,BG12#A) 

100  CONTINUE 

c ••*«•*£  have  COMPLETED  READING  sOUNp,*;  aNq  C«TS 

PRINT  1*>0 

l6o  fORpAT  ( 1H0,2*HSH$  CARO(S)  RE  PROoUCED  — ♦ / > 
oEao  1E5* (0(1) *I»l*NMKUaS) 

1 U5  rOPMAT  <8f10,6) 


pRInT  1T0* (B( I) * I*1*NMH0»S) 

1 To  fOR*mT(1H  ,8G12.M 

SET  NROw * 8 ( , 1 *ISTYR£(.l 

nHO#*N*RO*S 

pO  <f~0n  JJaltNUMVAR 
,|l*KLO(  JJ) 
i2skR0  < JJ) 

TF  ( Ji  ,EQ.J2>  UC  TO  VooO 
KlROWaNWOB*  1 
<}000  cUNTTM'E 

1 1 snmROMS* 1 
no  R:l'1  I*i  1 *nwOf 
R(I)fl  . 
tSTypE(I)s-1 

9010  fONTINUE 

rtDD  SLACKS  TO  COEFFICIENT  MATRIX 

m£LEm*" 

MCoLafl 

nO  9 T 0 o I * 1 * NRo* 
n£LEm«\ELEm»i 

ncol*ncol*i 

T A ( NtLfM ) » I 

A (NELEH)al. 

I A(nCOL**NEU£M 
9 1 0 0 CONTINUE 

I A(NCOL»l)aNELEM*l 


FILL  IN  COEFHCIENI  MATRIX 


» 
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mCOGUHsI 

nO  <**01  JJ»l»NuMVAR 

jl*«LO( jj) 

, *5  _ O r>  . i ik 


jl*«L0( JJ> 
j2*kRd< JJ) 

IF  (Jl.LT.J2) 


OC  TO  V300 


pEAO  92 5o*  (YTEMP(I) ,l*i *NMROwS) 

9250  fORmAT (8plo.6> 


9270 


rO  927*  I*l«NMNUhS 

yTEmP1»YTEmP<I> 

tf  <ab«<yTemRi) .le.ztolzE) 

nELEm*nELEM*1 
T* (NcLFM)*1 
A <NElE" ) -TTEMP' 

CONTINUE 
nC0|_*NC0LM 
lA(nL0l*1>*NELEM*1 
r,0  TO  9 A 0 0 


GO  t<3  9?7C 


9300  no  939*  J*J1*J2 


93r 


93p0 


9390 

«A00 


500 

25c 


575 


285 

?80 


9200 


760 

361 


nO  9300  I*i,nmhO*S 
CALL  GETPHI (I,jJ,CUTS«J) *ATEMP) 

IF  (ABS  (ATEMP) .LE ,ZTOLZF)  GO  To 
nELEM*nELEmM 
lA(NrLEM)«I 
A (NF|_EM)*ATEMP 
CONTINUE 
nELEM*\ELEM*1 
Ta (NfLFM) *NMHOwS*NCOG^8 
A <NFi  Em)*1, 
nCOLsNCOL*1 
la (NO0L*1 > *N£Lem*  1 
CONTINUE 
nCOGUB*NcOgUB*i 
CONTINUE 

tF  (NrLEM.GT.MAxA)CALL  ERR  (3) 
pEAO  2o0, <VARN*M( I) ,I*1,NUMVAR) 

FORMAT (16A5) 

pEAO  25 0 , ( PROBnA ( I j ,1*1,8) 

FORMAT ( 8 A 1 0 > 

*****[)ONE  READING  IN  DATA 
PRINT  575 

formaT(iho,7imfor  yoijh  information 
I ki 0 kRo,/) 

OC  28oi»l,NUMVAR 

P9INT  285,1,KL0(I) ,KR0(I) 

FCRmAT(1M  »3.5x,Ib»19X,I5,8x,l5) 

CONTINUE 

••••••SET  UP  STARTING  BASIS 

Oc  92q0  Ial.NROW 
JH(I)«I 
K IN0AS ( I ) ■ I 
tF(k5.E0.1)CAll  SCAR 
7F (k2,kE.O)360,JtO 
PRINT  36l 

fORmaT(1H0,5SHPAC«EO  MATRIX  BY  COLUM'S, 


variable  NUMBER 


ROW  NUMBER  below  Each  elem 


a 


■>66 

367 

365 

370 


?90 


30C 


r 


980 


<S9o 
1 0 0 o 

1 0o5 


1008 
1 0 1 0 
1 0 1 £> 


ran  »«  k «*s»  jmotwWB 

mOg  oof>Y  FUflMLSHSP  10  SBfi  ~ 


ifNT) 

tai.nelem/iui 

pO  365  1 s 1 . 1 X I 
kK« ( t-1) *1 1*1 

tF (kk.gt.nelehujc  to  3fo 

k*kk*i  -i 

tF (k.Gt.nELEM)r*kELEM 
PRINT  366»  (A(  J)  *vj*KK.K) 
p«InT  367*  ( I A <j>  ,j*kk*M 
FORM « T (1hO,11Gi2,4) 

F0om*T(1M  ,1111?) 

CONTINUE 
rONTINuE 
PRINT  ? 90 

fORw4T(1h1,19hS1aHTING  TO  ITERATE) 

PRINT  300 

fORmat  (1M0, AX,  iJT-STAGt.P^OttLEM.TX.j^LO-E1?  POUND, 9x,  llKJPPER  BOUND 
1 ,6X,i8^BRAnCH1nG  VARIABLE*//) 

•••••ready  to  s f ART  LOOP 
•••••ready  t o sJaRt  loop 
•••••re ad y to  sIaRt  loop 
•••••rfady  to  sURt  loop 
Dr  98o  I*1*LSTmAX 
I uRVR ( I ) *0 
Zl  slh ( i ) * i «e  so 
p ARfn  T * 0 . 0 
1 STkR ( > )«0 
LSTkl«D»0 
STgpRBsO.O 
1 ST  nUm»0 
nLPs  * 
tBRpar*o 
K.iOLFT*  1 

DC  9901*1 ,Nl MVAR 
KL  ( I>»MO<  1) 

K« ( I) «xRO  , 1 j 
PRINT  i005*STGPRh 
fORmaT ( Iho, 6X  *f?. 1 > 


mOi  ft»nolft-i 

NLP,  ' " 


3 K'L  P ♦ 1 

IF  (NLP. GT.maXLP)  call  LHRU> 


.iKSKRO  (NUMVaR) 
pO  63o?  J3 1 , JK 


t>300 


V’  ( J ) * o , 0 
CALL  LINPRg 

TF ( K 1 .fQ.1 ) CALL  lPPR(«1 
tF(lflg.ne.i»go  TO  1 0 1 v 

lF(NLP.E0.1)CALL  ERR<S) 
print  1008 

FORM/iT  (1M*,2AX,1imLP,  InFEaS.) 

GO  TO  5000 

T F (VAL-0U8TOL.Lt, BUB) 00  TO  1020 
PRINT  1015 

FORMAT (1H*,25X,VkLB  GT  BOB) 

G 0 TO  5000 

•••••PUT  This  PRC8LER  On  the  LIST,  ftnD  the  next  Empty  spot 
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j^pt^jjypa^HSDioaee ^ 


io?o 

102* 


10^0 


I0S0 

1030 

C 


lc*»0 

c 

c 


70u5 

7010 

701& 


r F (LSTNUM.GT.lSTmAx)CAlL  ERR<6) 

P«Int  |025,VAL 
fQR**AT  ( lH*t25»*&i2.4> 
i STnUM«lSTnUM*  1 
?LS  T» O (LSTnum) .iTGPRa 
7LSTPA (LSTnum) »HtR£NT 

jf  ( iaBC4fi.N£.c  ) LSTkl  (LSTnum)  «KL  ( IBRo  •») 
tFiipROaR.nE.O) lSTkR (LSTnum> »*R ( IbRo -B) 

7LSLB  <LSTNu«l  *vAl 
•••**nCT  OONE  «1Th  LIS'  rET 

•**«*C^E*tE  *CCU£0*  CU^ACTIFIEo  vgosION  OF  *— ELEMENTS  OF  XrO^ED 
•••••A-E  X,  IF  A CORBEPPONOS  TO  A LI-FA*  Va«t AdLE«  0«  ELSE  T-E 
• of  Guofc-EO  VARIABLES 

DC  1030  I«1,MJmv»h 
IF <KL0( I ) .NE.KBO(i) ) GO  TO  t0A0 
I!«KLO<I> 

XCODEO ( I ) *• ( 1 1 ) 

GO  TO  103'? 

xoooEod)»o.o 


oo  iost  I I«I- . iz 

XCOOEU ( 1 1 «xLUUEO  <1 1 (I  I) *-U  TS ( 1 1 1 
CONTINUE 

•••••XCOOEo  C«£ATE0*  C*ECK  IT 
OC  1060  I»l.MJMv*K 
IfjKLO(I)  ,E(-.kPC(I)  )G0  TO  1 060 
IL*KLOd  > 

[R»KRO< I) 

X*L«CuTS(IL)-CUTTul 
X*R»CuTS  dM)  *CUT  fOL 
tF{XCOCED<D  .LT.aXUCALL  ERR<7> 
tF(xCOCEO<I>*OT.>XR)CALL  err ( 7 J 
(-OntTNlE 

•••••now  «E  "lit-  tRT  10  FIno  an  UPPp;  dOuNO  »NO  ALSO  a 

•••••gSANCnlNG  vaRjaslE 

tFeaS.- 

rlFFMx»o.O 

TtJPVP  (lSTNum)  ■* 
mB*1 ,E70 

OC  2000  IRO*«l ,nmR0«S 
lFdCHK(IAU»).EQ.l)GO  TO  2n00 
RrKVAL»0. 0 

00  2'lC  I ■ 1 »NUM VAR 
•••••SET  InOEX  AnO  FRaC 
T»«kLO( I) 
tZ*kRO(I) 

tF<Iw.fq.IZ)700S,7010 

TN0E*«-1 

FrAC«0.0 

r, 0 TO  7900 

TM*I7"1 

nO  7*1«  I JsJHi  in 
TR»I,i*l 

TF(cLTSdJ),LE.AcOOEO(l)»AND.CUTS(lK,,GT.XCwnEO(I)>GO  TO  7025 
CONTINUE 

xXX«;UTS(Id)-XcOcEod> 

xAX«ABf(XXx> 


I 
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H 


i, 
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» 


TF(XXX,LT.CUTT0U GO  to  70 16 
XXX.XCOOEO(I)-COTS(IZ) 
xXXbAB? (XXX) 

iF (XXX.LT.CUTTOL)GO  TO  7017 
C*LL  E»R<7) 

70?5  pHAC»(xCODEO(I)-COTS(10) ) / (COTS ( IK ) -r UTS ( I J) > 
tN0Ev»I J 
r 0 Tl*  7900 
701?  tNDEX*TW 

F^ACbO  . 0 
r0  TO  7900 
7017  IN0FX«IZ 
fHAC,0.0 
7900  CONTINUE 

OlFFaO.O 

DC  2020  IJ»1W»IZ 
I I»I j*NKO« 

INOCOL*OA(II) 

lNDNXT*l»A  ( 1 1 ♦ 1 > "1 

do  2030  IIl.lNOCr'L’lNDNAT 

IF  1 1 A ( i j i ) ?nE.iroW>(>0  Tu  2030 

IF (INOEx.NE,-l)nTFF»DlFr.W(IJ)#A(lII) 

IF (IndEX.E&,IJ>DTFF»dIFf*< l.-FRAC)  • A(Ul) 

IF  ( INDEX. EG.  - 1 > Rr>WVAU»RuwVAL*XcODED  ( I > *A  ( I H ) 
IF (INDEX.EQ.I J)RoWVAL»RuwVAL*(i.-FRAC)*A(III) 
IF ( ( 1NdEX*1) .Efl.T J)DlFF»nIFF*FRAC»A( III) 

IF  ( (IN0EX*l)  .EQ.t  J>  #0SV*i  >ROwVAL.4FrAC^A(iii) 
2030  CONTINUE 

2020  COMlNOt 

rF ( ISTrPE  < IRO*)  .E0.-1 >D1FF«ABS  (qIFF) 

IF (0IFF.LT.U1FFMA*UIFfT0)6o  TO  20lf> 

DIFFmAbDIFF 
fLag(LSTNUm) "INDEX 
FLAfi(LSTNUM)»FLAC-(LSTNOM)*FRAC 
IBRVR (lSTNUM) b1 
2 0 1 0 CONTINUE 

IF (ISTYPE(IROW) ) 2040 12050 »2060 

2050  UP«ROwVA(. 

GO  TO  2000 

2060  TF  (B  ( IROW)  ) 2061  *2062»20<>3 

2061  tF<ro"vAl.|_E.  <0(IHOW>*<1.-FEAST|_>  > )Go  TO  ?0un 
TFEA^al 

r- 0 TO  2000 

2062  iF (pOwvAL»LE.FeASTl) GO  TO  2000 
tFEASbI 

r0  To  2000 

2063  tF(PO«IVAL.LE»  <8  < lf>OW)  •(  1.*FEaSTl)  ) ) So  TO  20«o 
jFEAsbI 

0°  To  2000 

2040  IF (ABS <B ( IHC4) ) ,tU.y.0>GO  T°  2p70 

XXXbI ,-ABS (POBVAL^B ( IRO«> ) 

IF  (ABS(XXX)  .UT.FtASlDGO  Tq  200  < 

IFEAS.I 
GO  TO  2000 

2070  IF  (ABS(ROwVaL)  .LEiFEASTUGo  TO  2000 

IfEASbI 

2000  CONTINUE 

C •••••OC'NE  — *E  HAVE  PICKED  A BRANCHING  VaRIaoLE  AND  STORED  IT  ON 
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c •••••The  LIST  v,HILe  TESTING  FOh  FEidlBlLlTl 
tf<ifeas.ne.i>gO  To  30O0 

PRINT  ?005, jhWVR(LSTnUm) 

20  0^  fORmTT  (1h*.*9x,*hN0NF»16A,i6' 
rO  TO  5000 

3000  PRINT  3005.UB* IBRVR(LSINUM) 

5Q05  FORMiT  <lH*,45x,Gl  2.*,  12X»I6) 

IF(IXPRIn.Eq.1)CaLL  XP«InT (XCOOED) 

IF (UB.LT.BUB) GO  TO  30IO 
pO  TO  5000 
30  1 0 r08*iiB 

DO  3020  I«1,(NUmVAH 

c *****now  begin  branching  procedure 

3020  XREST(I)»X<.C0ED(I> 

<5000  tF(nOLFT.EO.O>GO  TO  50B0 
c ***««S0LVE  NEXT  PROBLtM  IN  THIS  STAGF 

kL  < IoRPAR) »k8L  t NCLfT ) 
kH(I8HPAR)»KBR(NCLfT) 

<;TGPRB»STGPRB*.l 

rO  TO  1000 

c ****#WE  ARRIVE  HERE  if  "E  ARE  DONE  WtTH  A SI*gE 
5050  NUM»* 

rLBcI »E70 

DO  5060  I »1 « LSTNUM 
IF (ZLSLB(l) . GE.BLB) GO  TO  5n60 
NI;M«  I 

8LB»ZLSLB(i) 

5060  CONTINUE 

tF(hLB.GE.1.E7q)CALL  EkR(8) 
tF(hLh.GE.BUB-dOnTol>UU  TO  8000 

c *#**«NUM  is  The  entry  on  The  LIsT  ON  WHICH  «F  ARE  TO  BRANCH 
5063  fORmYT(1H0,20HdUnE  WjTH  THIS  STaGE) 

TF  (IbRvR  ,NE.O)GO  TO  506A 

CALL  ERR (9) 

506*  PRINT  5o63 

PRINT  R065tBLB,  BOB, 21.51  NO  (NUM)  ,iBRVR ,NUM) 

S06b  FORMAT ( 1 H ,6HBLBe  ,GiZ.4*8H,  BLlB»  ,Gi2,4*£?H*  BRANCHING  ON  PROBL 

ipm,fc. 1* i7h*  variable  number. i6) 

TF(kA.fO,1)PRInT  50651 

6n651  FORMAT (1H0,32h«*********PR£SENT  STATES  Of  LIST ) 
tF<K4.FQ.1)PRInT  50652 

5*652  fORmaT (lH0.6HPRObN0.5X*6HPARENT,5x,XuLlSTKL»SX.6HLlSTKR.5Xf 
1 1 IHLOWER  BOUND, UX.12HBHANCH.  VA«. ,CX . 4HFLAG*/) 

IF (K4.EQ.1 ) 50653, 5065? 

6*653  nO  5?oo9  IbI.LsTNUM 

print  c0001 .ZLsTnOII) .ZlSTPA u) , LSTkL»1>»  lsTKR(I) 

60OO  I fORmaT  ( 1 H ,F6.l,5X,F6.1,i»X,l5,6x,I5) 

(ZLSL9 ( I) .GE. I.Eto) 5QCU2, 50004 
5«cr2  pRJnt  50003, IBhVR (II ,FLAG(i) 

5 -cn 3 rOBwAT (lH*,48X,3ROFF,4X,9X,l4,9x,Flo.3) 

<*0  TO  50004 

p.3i*  6"[nt  *ooo5,zlslh(d , ibrvr < i ) , Flag ( 1 1 
« , b« * T 1 1h*,**a,G1 1.^,Vx,I4,9x,Fi0.3) 

* -OnTI*,'  I 

••!«.»  «o«66 

• * r *•!* ' !*••» 

•« * •»,,»  c 
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TBRPAH»I8RV«<NuH) 

••***now  we  know  the  parent  and  branching  Variable  for  the  nfxt 
•••••STAGE— SET  IP  2 UK  3 NE*  PrQ8LF«S 
•••**F III  KL  AND  KR  VECTORS 
DC  5070  I«1*MJmVaR 
KR(I)»KRO(I> 

KL (I>«KLO(I) 

7N8ACK»ZLSTN0(nDk) 

DC  5080  I » I * NUm 
I I»NUM“I ♦ 1 

If  UnBACK,E.G.0#0)6o  to  si  1 0 
lF<ZLSTNO(llI ,NE.ZNBACK)QO  TO  5, -80 
DO  5oV('  Ik«I*NUM 

IF(ZlSTNO(IN).EQ.ZlSTpA(jp )GO  To  *o95 

continue 

III»IBRVR(1k> 

JF (.not. ( (LSTkL(II) .6E.KL  < 1 1 1 > ) .AND. (L*TKR ( 1 1 > .LE.KR « 1 1 1 ) ) > ) 
1 GO  TO  5100 

KL  (III)wLSIkL(II) 

KR  ( 1 1 1 ) BLS  f KR ( 1 1 ) 

ZNBACk*ZLsIfA(II) 

continue 

rONTINLE 

•••••NCW  H4Ve  TU  DIVIDE  UP  The  K-SET  FOR  THt  rRaNChING  VaRIarLE 
**««*hut  first,  remove  The  parent  problem  FnoM  The  list 
7LSLB (NUM)»1.E70 

*****SfT  Up  Two  OR  THHEE  PROBLEMS 

if  ( <Cl»G<NUM)  .L  ! .KL  dBHPAR)  ) .OR.  (FLA^(nuM*  .v>T.kH<  IBRPAr)  > ) 

IfALL  EQR  < 1 1 ) 

•••••CREC*  to  SEE  IF  FLAVJ  PRtClSELY  rOUAL^  ROME  COT 
I »*KL ( IBRpAs) 

IZ*kR( IBRPAR) 

DO  5120  J«iw#IZ 
Z JW  J 

XXX»Zj-FLAU(NUM) 

XXX»ABS(XXX) 

IF (XXx.L£.CUTTOL>00  TO  5l3n 

continue 

Ta*FLAG (NUm) 

tF(ix.fq.kl(ibrhar))go  td  51*0 

KbL  < 1 > »KL  < iBRPAH) 
kBR ( 1 ) a IX 

K^L  (2) »ix 
kBR(2)«ix*1 

IF  < ( IX* 1 > .EQ.KR( IBRpaR) >5150.5160 

K'0Lft»2 

G 0 Tr.  *000 

KBL  (•»  > » I X ♦ 1 

kBR ( 3 ) «KR ( IBRPAH ) 

mOLFTw'’ 

GO  To  6000 
KBL  (6 ) ■ I X 
kBR ( 1 ) »I X*1 

kBl (?) «ix*i 

kBR  (2) «KR  ( IBRPAH) 

nOLFT»2 

gO  To  6000 

IF< (J.EQ.KL(IBRHaR) > .oh, (J.ED.Kr(IBPoaR>  > >C«l  L EHH,l0) 


r 1 
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K«L(l)»KL(IbRPAK) 

KbR(l)«J 

k«L(2)»j 

KBR(?)«KR(lbRPAR) 

nOlft«? 

6000  iXX.STr.PRB 

STGPRB-IXX 

sTgpRh*stgphb*i • 

rO  Tr  5000 

C **#««OONE“-  PRINT  OUT  TrtE  RESULTS 

9000  pRIn!  9010, (PRO^N4 ( I > *I«1 *8) 

8010  F°RmaT ( 1 H 1 8 A 1 o ) 

PRINT  0020, BUS 

8020  fORp;T (lHO,31HOBJtCTIVt  FUNCTION  AT  rPTImUM  ,G12**) 
PRINT  H030 

8030  FORMAT (1H0,2BHVAFIABLE  VALUES  AT  OPTtMUM— ) 
rALL  XPRINT (XBEST) 

FNO 


SUBROUTINE  XPRJNT (Z) 

cOMMON/FlRST/KLOdOO)  *kR0(100>  *KL(ln  > *KR  ( luo ) »XCOotO  1 1 00  > *X0EST 
1 (100)  , mill  00)  .CUTS  ulOO)  »ZLSTNO,7oO)  . ZLST?A  i toO  > ",  LSTKL(700)  , 

? LSTkR (7oO ) »ZLSLh (7oO> ,IbRVR<70o) *FHG  < 7og>  *K8L  <3>  *KBR<3>  * 
3vARn«M ( 100) , PHOBNA ( 8 ) tMA*VAR*MAXCUT,i  STRaX,m»XRO*#mAXA, 

4 NMRO*5»NUmVAH,ICHK(1O0>  *VAL*LFlG 
oIMERSiON  2(1) 

PRINT  )0 
10  fORmaT(IHO) 

I L*  1 

20  tUsNUMVAR 

jF(  (nUmVaR-IL)  .(jT,7)  IU*IL*7 
PRINT  40, ( vARnam ( I ) , i»1L» IU) 

40  fORmAT (1H0,3X,aS,7 ( 1 2A, Ab) ) 

tz«iu-tl*i 

PRINT  45 • ( Z ( I ) ,I*lLtIU) 

45  fORmTT(1h  ,Gl2.4,7 (SXtG1*.4) ) 

TL»IL40 

tF(IL,LE.NumVAR)gO  to  20 
PRINT  10 

return 


SUBROUTINE  LPPR(Y) 

COMm0N/FIRST/klU(100>  iKRO(IOO)  ,KL(lr  '•)  fNR(lwo)  * XcUptO  (100)  *X0EST 
) (ion) ,4(1100) ,CUTS( 1100) • ZLSTNO ( 7oo ) .ZLSTPA17  00) i lSTkL(^OO) , 

? LSTKR(700) ,ZLSLb (700) • I«RVR(70o) *FL«G(?o6) #k0L(3) ,KBR (3) t 
7VARN AM (100) »PRObNA (0) »maXVAR»MAxCuT»i  ST MAX*n»XRO«t *MAXA , 

4 NmRORS*nUmVAR,IcHK(10O) *VAL*LFlG 
nlMENSTON  Y ( 1 ) 

PRINT  10 

FORMTT (1H0,29MPAC«ED  UR  SOLUTION*  I X(I)> 

I4«KRo<NUmVAR) 

00  30  I>  1 • Am 

T F ( aBS ( Y ( I ) ) .GE.l ,E-10) PRINT  20,1 *Y (I ) 

FoPMATdH  *l5x,l6*2x,Gl0.4, 

continue 

PRINT  40 
FORMAT (1M0) 

RETURN 

FNO 
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SUBROUTINE  ERH  ( I > 
pRInT^O 

format <*iProgham  mogg  aborted  because. • • .*> 

ROTO  ( 1"  It  102*  10J*104»I05»  106*107*  108  *lO?»ilMein>  *1 
PRINT  ?01 
C*LL  EXIT 
print  202 

C*LL  EXIT 
PRINT  203 
CALL  EXIT 
PRINT  204 
CALL  EXIT 
PRINT  205 
CALL  EXIT 
PRINT  206 
call  EXIT 
PRINT  207 
call  EXIT 
PRINT  208 

call  exit 
print  209 
CALL  EXIT 
PRINT  210 
call  EXIT 

PRINT  211 
CALL  EXIT 
RETURN 

FORMftT <*oVaRIAbLE  CARO*  OUT  OF  ORDER— LOOK  "FA*  MOqG  LaBEL  105*> 
format (*0MAXCUTS  EXCEEUEO— look  near  HOOg  L*»mEL  117  OR  124*) 
formaT(«omatrix  a exceeded— loo*  near  rOqg  label  ?aoo*) 
format <*oLPMAX  EXCEEDED— LOOK  N£AR  MoGG  LABEL  1005*1 
fORmAT (*0lNlTlAL  LP  INFEASIBLE--LOOK  near  MWGG  label  1 00®* ) 
FORMAT(*oLlST  LEhGTH  txCtEDEU— LOOK  mEAR  MOwr  LABEL  1020*1 
FORMAT (*0XcODEO  VIOLATES  CUTS— LOOK  nEAR  MOwr  LABEl  1C60  0r  7025*1 
FORMAT (*ONO  BRANCHING  NODE  FOUNo— LonK  NEAR  mOGG  LAbEL  5060*1 
FORMAT, *0NO  FEASIBLE  POINT  FOuN0..LOnK  NEAR  mOGG  i-ABEL  5064*, 

format i*qNo  bhanching  possible  on  Variable  lmoSen--look  near  mogg 

1 L^BFL  5lj0#) 

FORMAT  <*OFlAG  computed  impropeRly— L^OK  nEAk  mOGG  label  5Ho#> 
fno 
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SUBROUTINE  SC A i L 

rOHMON/FI«ST/KLO(loO) tK«U UOO ) • *L ( 1 0' ) tKR < l«n ) • XC00ed  ( 1 00 ) tXgEST 
1(1 00) • « ( 1 100) (CUTS (1100) *2LSTNO(700) .ZLfTPAlTOO) * LSTKLC700) , 

? LSTk«(700)  ,ZLSL8  (700J  tlBRVR(TOo)  .FL^ffoO)  »KBL(3)  .Kapil) t 
TVARnJm (100) tPHOBNA(B) ,MA*VAR,MAxCUT ,t  STMaX ,MaXRO*  ,MAXA , 

4 NMROtaS  t NUmVAR,  ICHk  ( 1 00  ) * VALtLFLG 
cOmmON/WORkI/  B<35o) ,A  (3So) »Y (3S0) ,YTFMP<35ut ,A (5000) »E  C 5000 ) * 

1 I A ( 5r  00  > , IE (5000 ) ,LA<iJ02) . LE ( 200?) . I?NAM t , 3o2 ,2) ,KInBA§ ( 13o2) f 
? JH<350>  tlSTYP£(350) iNAME(20>  *nTEMP<20>  *CM1n,C0N0,ERHAX* IFFEZ, 

3 lNVFRQtlOBJ»lHO»P.I rCHilTCHA*lTCNT.lT^F^CI»TVIN»IVOUTt JCOLP.KINP, 

4 XST AT, N«OW.NcUl* NElEH.nET A *NLELEM,mlET A tN»FLEM*N6ETAtNUELtMt 
* NUET4»SUMINF,K3 

f-OMMON/ BLOCKS  ZlOLZEt2TOLPV»ZTCoST,NUMAX,NTH*X*NEMAX,aROtifMA,QBA» 

1 QF? . GEO*QBL»«PL*flMI iUA»QB,(JC*qE,0F.06»QM»«t  *UL*Om»QN.QO*QR,OU, OZ 
nO  ino  IXX»2,NmR0«S 
RHAlL»1 ,E7o 
bIG»-1.E70 
last«la (NC0L*1)-1 
iFIRST«LA(nR0W*1) 
nO  ?00  IXY.IFIRSTtLAST 
tF  < i a ( t xy) . ne»i*x)goto2oo 
IF ( A6S (A ( IXY> ) .LT,SMALL)SmaLL«AbS(A(TXY) ) 
tE(ABS(A(IxY)) ,UT.BIG)bIO»ABS(A(IXY) ) 

200  rONTTNUE 

aV.SoRT (SMaLL*0IG> 

7L?Av«AlOG(AV) /Al0g(2») 

L2Av.^T(ZL2AV) 

TF (ZL2AV.LT.Q..AN0,L2AV-ZL2AV,Ge.,5)I pAV«L2Mv-1 
IF(ZL2«V.GT.0..AND.ZL2AV-L2AV.Ge..5)i 2av«l2mw*i 
nlV«^.**L2AV 

OO  300  IXY.IFIRSTtLAST 
TF(IA(IXY> ,ne,iax>gotojoo 

4<1Xy)»A<IXY)/d1V 
300  CONTINUE 

R(IXX)«B(IxX)/oIV 
100  CONTINUE 
retupn 

FnD 
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SUBROUTINE  LlNPNG 


3VARNAM (100) »PR0BNA(8) *M*AVAR*MAxCUT.i  SIRA* »«aXROw »mAXA « 

4 NmROWS  *NUmVAR, ICHk ( 100) i VALiLFlG 
common/ WORK  1/  H'»35Q>  » A (3*0)  *Y  (350)  ,ytfMP(35u)  , A(50o0)  «E  (5000)  • 

1 I A (5oOO> ,IE(50C0) ,LA(1J02) .LE(2oo?t . ICnAH i, 302 ,2) .KInBaS ( \ 302) , 

? JH(350>  *ISTYP£(350> *NAME(20) »nTEMP(20>  »CMIn,COnO,ERMAX. IFFEZ, 

7 lNVFPQ,IoBJ.lRO«P.lTUH»ITCHAtiTCNT.lTRFRQ*TVlNtivOUT,jCOLP,KTNP. 

a xstat,nR0w,ncol.nelem,neTa,nleuem,mle!a.nwflem»ngeTa,nuelem, 

s NUETa .SUmINF,K3 

rOMMON/BLOCK/  ZTcLZE*^1 OLPVt2TCoST,NoMAA,NTMAXtNEHAA»«RO»WMA.QBA* 

1 OFT,'.EO,QBL»0RL*QMI,UA,Qb,Uc,(JE,QF.OG»QM#'-TtOL,«MtQN,OO,QR,QutOZ 


T TcnT«o 

tTch;*" 

Sr T UP  STARTUG  BASIS 
nO  R 1 0 0 J*1*NC0L 
kINHaS(J)*0 
nO  roo  I*1,nR°4 
TCOL-vlH  (I) 
mLsnROw 

OC  300  K* 1 f KUmVAH 
MR3KL <K> ♦nKG* 

IF( (ICOL.gI.RL) .and. (ICUL.lT.MR, )G0  TO  700 
ML»KR (K) ♦nKOW 

continue 

fFUCOL.GT.MDGU  To  700 
r, 0 To  <300 

JH  < I ) s I 
CONTINUE 


1000  CALL  INVERT 
tTSINV  * 0 
CALL  ITEROP(O) 

: Simplex  c*cle 

”l 500  CALL  FORMC 

CALL  Sh1FTr(3»4) 

I TCh*0 

1 TOO  CALL  BTRAN 
CALL  PRICE 

tF  (CMIN  .LE.-ZlrOST)  GO  TO  3oOo 
TF  (XSTAT  ,EQ.  «I  ) GO  TO  200o 
XSTAT  » QBl 
fiO  TO  6000 
2000  mSTaT  ■ ON 
r,0  TO  6000 

3000  cALL  UNPACK (JCOLP) 

CALL  FTRAN(l) 

FKMAX*r'  • 


013 


riO  8000  I"1»NR0* 
pkM4X»ERMAX*Y ( I ) «YTEMH ( I ) 

8000  CONTINUE 

nlFxX.CMIN.ERMAA 

niFxx«#as<oiFxx) 

TF  (nlFXX.lE.ZTCoST)  00  TO  850o 
fF  (k3.NE.1>  GO  TO  6^0 
pHlNT  9500,CMIn»eRHAX 

95oO  fORMAKIH  1 1 OX t 6hChIN“  tF  16,8,5<#7hfoMAX«  *M6.8) 

8ioo  if  (ERnax.ce.o.)  go  to  asoo 

IF  (ITCH.GT.O)  go  TO  1000 
tTCH«JC0LP 
tTChA«iTCHA*1 
CaLL  SHIFTR (A, 3) 

GO  to  1700 
8500  CONTINUE 

CALL  CHU2R 

tF(XSTAT.EQ.QU)Go  To  8000 
TVOUT»JH(IROWP) 

tvin  ■ jcolp 

call  UpBETa 
kIN8AS( JCOlP)  ■ iBOWP 
kIN8AS(IV0UT>  ■ C 
JH(IROWP)  » IVIN 
tTCNT  > ITCNT  « 1 
tTSImV  » ITSINV  ♦ l 
CALL  I TEROP ( 1 ) 

IF  (NElEH  .GT.  8000)  GO  TO  1«00 
CALL  WRETA 

if  (ITSINV  .GE.  INVfRUJ  go  TO  1 qQO 
IF  (ITCNT  ,GE.  ITRfRQJ  GO  To  6®0O 
GO  TO  ’500 

6000  CALL  I TEROP ( 1 ) 

SET  PaRmS 

nO  7noc  I“i»nro* 
jHX»jh<I> 

IF  (JHX.LE.NROw)  GO  TO  6500 
W<JuX-NROw)«X(I) 

8500  CONTINUE 
7000  CONTINUE 


VAL»-X ( IOBJ) 
lPLG*1 

IF  (vSTAT.EQ*QBL)  LFLGaU 
P«IN]  9600, IT CHA 

96o0  fOPM»T(IH*,108x*i8hSTABILITY  COijNT  b tl^l 
7100  continue 
return 
fNO 
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PBOii  QOiiY  EUKHISHEU  TO  DDC  - ■ — 


SUBROUTINE  FOMmC 

c 

rOMMON/WORKl/  B O50) .A(3!>0) t Y (350 ) »YTEHP(35uj ,A(S000> *E(50Q0) * 

1 IA(5nOO) , IE (5000) ,LA(]J02) .lE<2002> » ICnAM l , 302*2) .KINBA$  < n02> , 

? JH(350) tISTYPt(350) *NAmE(2o) »nTEmp(20) ,Cman,COND,ERmaX»IFFeZ, 

•»  INVFaQ.IOBJ»lHO»P*lTCH*lTCHA»iTCNT. iTRf^Q»TVIN.IvOUT. JCOLP.KjNP. 
4 XSTAT»NROtt»NcOL*NELEHfNETA.NLELEMtNLE!A*Nv»FLEH,NGETA.NUfcCEM. 
p NUETA«SUMINF,K3 

c0mm0N/8L°CK/  2lC’LZE.^[OLPv,2TCOSTtNHMAA,NT«AX*NEMAA,QRO,t!MA,OBA, 

) QFI .GEO»QbL*QPL *QMI»WA*ab.UC»QE#OF.OG»OM.wI,ULtUM*UN*QO*aR,OU,QZ 

c 

c 

xST  AT*CF 
iFFFJ  - 1 
nO  1*0  I « ItNROta 
V < I J « 0. 
loo  rONTlNLE 
SUM  ■ 0» 

c 

r 0 loon  I « lfNHO* 
tCOL  * JH ( I ) 

If  (ICCL  .Gl.  GU  TU  500 

IF  (ISTYPE(ICOl>)  200  * 1 00 0 * 5oO 
C 

2oO  tE  ( ABS(X(I))  .UE.  ZTOL^E)  60  to  lo-O 

tF(x(d  .lt.  o.^  Y<n  x *i. 

tF  < X < I ) ,GT.  0,)  Y ( I ) * -1, 

SUM  s SUM  ♦ Aa$(X<|) J 

r.O  TU  510 

c 

500  TF(X«I)  .GT.  -ZToLZe)  GO  TO  loOo 
Y<I)  ■ 

sUm  * SUM  - X(|) 

510  TFFF7  ■ 0 
XSTat  ■ QI 
1000  rONTTNt'E 
C 

sUminf  * sum 

IF  I'rFFEZ  ,LE.  0>  go  IU  *000 
Y<IOnJ)  ■ 1. 

C 

9000  RETURN 

fno 


THIS  PAGE  IS  BEST  QUALITY  PRACTICAL 
yiflOil  OOPY  FUBMSHH)  TO  DDC ^ 


SUBROUTINE  BTRaN 
C 

tOmmON/WORkI/  H 4350) » A (3*0) t Y <350 ) . YTEMP (3Su| «A(5000>  »E<S000>  * 

1 IA<5n00) , IE (5000) »lA U 302> »LE(200?) .ICnAhu 302*2) *K1nBA$(130?) , 

? JH(350) *ISIYPE(350) *NAHE(20) .mTemp(20>  ,CMiN#COND,ERMAX.lFFFZ, 

•*  lNVFR0*l08J»lH0«P*IICH*ITCHA*lTCNT.lT^rftQ»rVlN»IV0UT»jC0LP,KlNP, 
4 XSTAT.NRoP*NcOL*NELEH,NETA*NLELEM,MLElA*NtopLEW*NGETA.NU£LrM, 

* NUPT4*SUMINF,N3 

rOMMON/BLOCK/  2lOLZE,ATOLPV,2TCoSTtNWMAA,NT«4X.NEMAX*ORO.«MA,OBA, 

T OF! »GEO»QBL*OPL*QMI*UA*QB*UC*UEtQr.QG*QH*wT ,«L*«M*ON.OO*QR,OU*OZ 

C 

TF  (NeTA  »LE.  go  TO  9000 

nO  1000  I • 1,NETA 

|K  > NFTA  - I ♦ 1 

Lt  a LE(IK) 

kK  a LE(I**1>  - ) 

TPIV  a IE<Ll> 
nP  • E(LL> 
nY  a Y(IPIV) 
nSUM  ■ 0, 

TF  (KK  .LE,  LL)  GO  TO  600 
lL  a LI  ♦ 1 
nO  5fl0  j a i.L*KK 
TR  a IF (J) 
nE  a E(J) 

nPROO  a OE  * Y < m > 
nSUM  a DSUm  ♦ oPROD 
5oO  rONTINUE 
C 

600  Y ( I P I V ) a <DY  - DSUM)  / OP 
1000  rONTINUE 
C 


9000  RETURN 

fnd 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
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* 

I 

♦ 


C 


c 


c 


SUBROUTINE  PRICt 

cOMMON/FlRST/KLOUoo)  *KR°<l00>  »KU<in;)  *KR(luo>  **C9oEO(100>  tXREST 
1 (100)  «M  (Uoo)  tC^TSdlOO)  » ZlSTNO  ( 7 o 0 ) .ZLSTPA1700)  , LSTKL<7*>0)  , 

? LSTKR ( 700 ) ,ZLsCb(7oO*  *IttRVR(700) ,FUG(?o6> »kBL (3) ,KBR(3) » 
3VARNAM(100>  *PROBNA(b)  *MA*VaR»MAxCUT,i  STMaA**axROiMmaXa« 

4 NmROWS*NUmVAH,ICHK ( I 00) *VaL»LPl6 
tOmmOn/WORkI/  H«350) , A(3?0) .Y(350) ,YTFMP(35«) ,A(5000) ,£(5000) , 

1 I A (5o00) , IE (500  0) ,LA(l302) »LE(2oo?i , ICnAM t , 302,2) ,KINBAS(13o2) , 

? JH<350> ,ISTYPE  <350> »nAME »2o>  #nTeMP »20> »CM4n,C0ND ,E«Ma* * iFFpZ , 

7 inVfRQ, IOBJ*IrC»<P,iTch* ITchA* iTcnT.iTRfRG,tVIN, IVOUT, jcOLP,KTNP, 

4 XSTAT,NRow,NcOL,NELEft»NETA,NLELEM,NLElA*NwFLEM,NeETA,NUfc.LEM, 

5 NUETa,SUmINF,N3 

rOMVON/BLOcK/  ZlOLZE,ATOLPV t2TC0ST,NPMAX,NTH4X,NEMAX,0R0,GMA,QHA, 

1 Qft»CEO»QBL,UPl*QMI,UA»QB.UC*QE,QF.QG*QR»Ut*OL.Qm«QN,QO»QR,OU,OZ 


jCOLP  ■ 0 
CMIN  ■ 1 »E 
pO  100o  j 


10 


l*NCOC 


jF(j  .LE.  NRO*  .AND.  15>TYPE(j)  ,NE.  i>  Go  Tw  1000 
tF  (KINBAS(J)  4NE.  0)  GO  TO  1000 
IF  (ITcH.EO.J)  go  TO  1000 


OL«NROw 

no  3:0  Ka  1 , NUMV  Ap 

qR*kl  U>  *nrop 

tF  (<J.GT»QL>  •AhO»  (J’LlfGR))  GO  To  1000 
qL»kR ( K ) ♦Nrob 
300  CONTINUE 


C 

tF  (J.gT.QU  GO  to  1000 
nsup  ■ 

lL  ■ LA(J) 

KK  a L ^ < J ♦ 1 > " * 
pO  500  I a lc*kK 
jR  a I 6 ( I ) 

PE  . A ( I ) 

pPROD  - DE  * Y ( I P ) 
pSUM  a OSUm  ♦ [)PP0d 
500  rONTINUE 

tF  (rSUM  .GE.  cMIN)  G°  T0  1000 
C*IN  ■ OSUm 
, „ JCOLP  a J 

1000  cONTtNUE 
return 
fnd 


Ct-17 


W 


tBIS  PAGE 

pRfiiiCOFX 


IS  BEST  QUALITY  PBA'CTlCABT® 
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SUBROUTINE  SHIfTr<IOlD*INEW) 

C 

rOHMON/WORKl/  0 050)  ,M3»0|  *Y(350>  ,YTEMP(35u,  , a (5000)  ,£(5000)  , 

1 IA(5n00>»IE<5000)  *LM1JQ2)  »L£<2002>  » ICn*M»  j 302*  2)  *KInBAS  ( 1 3 o2 ) « 

? JH(350>  *ISTYPE(35o) tNAME (20 ) *NTEMPf20>  ,CMAn,CON£),ERMAX  . lFFEZ, 

Y INVFRG* lOBJ, iROWP, I TCH t ITCHA*  jTCNT « I^FRQt  j VIN, IVOuT  * jCOLP ,KT  NP  » 
4 XS]aT  *NROIY»NcOL*NELEM,NETA*NLELEM,NjLETA»N«»FLEHtN(JETAfNUELEMt 

A NUETAfSUMlNF.Kj 

cOMHfiN/BLOcK/  ZjOLZE » 2I0LPV *ZTCoST ,NBMAX,NT*aX ,N£MaX» aRO»UMA,QBA* 

1 OFT  »CEO»QBL»QPL»QMI»QA»QB»UC*QEtQF.QO*QH»«T,gitQMtUNtQO»QR,QU»OZ 

C 

DIMENSION  BARHaY(U00> 

FOUIvAlENCE  (faARpAY(l) »B(1> > 

C 

IFO  . (lOuo  - l)  * NRMAx 
iEN  ■ (INEw-  1)  • NRMAA 
C 

dO  1000  I » I»NR0» 
rARPaYUFN  ♦ 1)  • 8ARHAYUF0  ♦ j) 
lOOO  c°NTINUE 
return 
FNO 

SUBROUTINE  UNPACK < IV) 

c 

rOMMON/WORKl/  B (j50) ,A|3»0) .Y<350) .vTEMPjSSu, ,A<5000) »E(POOO) , 

1 I A (5n00>  * IE  < 50c 0 > ,LA(1JQ2) iL£(200?' * ICnAM \ , 302*2) ,KINBAS<i3o?) 

? JH ( 350) » ISTYpE ( 350) »NAM£(20) *nTemp  f 201 » CM*n»C0n6, ERMAX 1 1 FFfZ* 

Y iNVFRO.lOBJf lR0«iP.IICH»lTCHAiiTCNT*IT^FRO»TViN*IVOuTtJC0UP,KINP* 
4 XSTAT,NROK,NcOLtNEu£M*^ETA,NLELEM,NLElAtN»»FLEM,NGETA,NUELEM, 

A NUETa*SUmINF«K3 

rOMMDN/BLOcK/  ZlOLZE*2TOcPV,ZTCoST,NBMAA,NTHAX*NEMAX»OROtUMA«Q8A* 

1 OFI*CEO*QBU*OPL*QMI*QA*ObtUC*OEtOF.OG»QH*wr,OL*UM*UN*QO*QR,OUfQZ 

C 

no  100  I ■ i«nrOw 

Y < I ) ■ 0. 
loo  r^TlNuE 

c 

LL  > LA ( I V) 

K*  • LMIV4I)  - 1 
nO  200  1 ■ LL*kK 
tR  ■ I * < I > 

VH»)  > A < I ) 

200  CONTINUE 

c 

RETURN 

pNo 


1 


Cr-18 


t 


SHIS  PAGE  IS  BEST  QUALITY  FRACTIONAL* 
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d 


SUBROUTINE  FTRaN(IPAR) 

cOmmON/wORkI/  R(35o) *X<3?0) »Y(3sO) ,YtEMP(35v>  , A (5000) »E(5000) * 

1 I A (5000) f IE ( 5 0 0 0 ) # L* 11302) *LE(200?1 t ICnAM \ ) 302 *2) «KI NBAS (130?) , 

? JH ( 350 ) • 1ST YPE ( 350 ) *NAM£(20) »nTEmP(20) »Cm1n*COND. ERmAX  » IFFEZ • 

•,  INVFROtIOBJ*lRO-P»lTCH*ITCHA*lTCNT«lTfiF^Q»TVINt 1 yOUT » JCOLP*  KTNP • 
4 XST AT »NRoX»NcOL» NELfcH.NET A tNLELEM.NLElAfNtoFLEHtNGETAfNUtLEM* 

1:1  NUETA,SUmINF,K3 

cOppON/BL°CK/  ZlCLZE»AT0LPV,ZTC0ST,NoMAX,NTM4X»NEHAX»QR0*UMAiQBAt 
1 OF  I »GEO'QbL'ORL  *QM I fUAtQB'OC'OE  *OF.OG*QH»uj *QL»0m*QN»00»QR.QU*QZ 

C,0  TO  (100,110)  tl^AR 
100  NFE  ■ 1 
NLE  ■ NETA 
go  to  ?oo 

no  nFe  ■ nleTa  ♦ i 

nLE  « kETA 

200  TF  (NFF  .GT*  N(.E)  go  to  vooo 
nO  lnoo  IK  ■ NFEfNLE 
LL  « LE<IK) 

KK  . LE ( IK*1 > - 1 

IPIV  ■ IE  ILL) 

OT  ■ Y(IPIV) 
pT  ■ 0Y/E<LL> 

Y I IP  I V ) ■ OY 

TP  <KK  «LE.  LL)  GO  TO  i00° 

LL  ».LL  ♦ 1 
nO  500  j « LL*KK 
TR  ■ IF ( J) 

Y(IR)  ■ Y(lP)  - E(J)  • OY 
500  continue 
1000  CONTINUE 
9000  CONTINUE 
RETURN 
FnD 
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SUBROUTINE  chuzr 

c 

cOmmON/WORk 1/  0(350) #A(350) »Y{350) .YTEMP(35u> , A < 50o 0 ) »E(S000) t 

1 14(5000) f IE (50O0) «LA()U02) »LE(2002t . ICn*M i j 302. 2) .K;n0A$ ( 1302)  . 

2 JH(350) .ISTYPE(350) .NAME (20) .nTEMP<20> .CMIn.COND.ERmAX.IFFFZ. 

3 lNVFR0.l08J.lR0*P.I)CH*ITCHA* I TCNT  » iTRpfto.j  VIN. I VOUT » JC0LP  «K1NP  » 

* XsTAT,NRo*.NcOL»NELEM,NETA.Nl.ELEM,NLE!A*N«»fl.EM,NQETA,NU£LEM, 

NUET4.SUmINF.k3 

tOmmon/blOck/  zTolze.xJolPv.ztcoSt.nqmax, nth^x. nemax. Qro. Rma* 004 , 

1 OF  I .GE0*0BL.QPL*flMI * WA *Q0 »0c *QE *QP .QG *QM.NJ .QL«GM»QN.OO*ttR .QU.QZ 

SELECT  Pivo!  ROW/VARIABLE  !o  LEAVE  IhE  OASIS 

7T0LCR«1.E-A 
ZTOLXXal .E-10 
xMinT»1 ,e!0 
XMIN2«1 .E^o 
XMin3«1 »e!0 
IROwP 1*0 

]ROWP2«0 
lROWp3aO 


nO  2000  I»1.NR0W 
JF  (ISTYPE(I) .E«.0)  GO  TO  2000 
tF  (ABS(Y(I) ) ,Lt.2T0LCP)  GO  TO  2000 
iCOlaJH (I) 


JF  ( (ICOL.LE.NROIx)  .AND;  (ISTYPE(l>  .LT.O)  ) 
XRAT IO“X ( I ) /V ( j) 
iF(XPATIO.lT.*ZIcLZe>GOTo2000 
TF(Y(I).LT,0.)G£)T02ooO 


tF  (JRaTIO.GT.XMINj)  go  TO  2000 
XMJnI axRAT  JO 
iROWplal 
r,0  TO  2000 


Gu  TO  I®00 


1000  IF  (A0S(X(I) ) .LT.ZTOLZE)  GO  TO  I50o 
xRATtO»X(I)/Y  C i ) 

IF  (XR4TlO.LT .0;)  GO  TO  2000 
I F (XR4TI0.GT.XMIN2)  GO  TO  200o 
XMlNPaxRATiO 

yR0Wp2aI 
r,0  TO  2000 

1500  xXXa4BS<Y<I>> 

xR4Tt0»ZT0lXX/xXX 

IF  (XR4TI0.GT.XMIN3)  GO  TO  2o0o 

XMlN3aXR4TlO 

TROwP^al 


2000  f-ONTINUE 


TEST  FOR  OUTGOING  VECTOR 

tro*p«iro*pi 

xRATfOaXMlNl 

TF  (XR4TI0.LE.XMIN2)  GO  TO  300o 
TR0UP«TB0IXP2 


C-20 


XRATIO.XmInZ 

c 

3000  tF  (XR«TI0,I.E.xMIN3)  Go  TO  400o 
TR0WP*lR0»lp3 
X«AT?0»XMIn3 
C 

*00?  if  (IROWP.LE.O)  XSTAT«OU 

t»irowp 

if  1«3 ,NE  • 1 > RETURN 

print  9000,IROwR,X(I) ,T (I) ,JH(I) 

9000  FORMAT  Oh  IROWP-  » !♦»  i!X,  6HX  ( I ) * ,F16  8*2X*6my<I>* 

] .7HJH(IJ«  ,I4| 

return 

fND 


»F16.8»2X 


subroutine  UPBETa 


c 


c 


i 

? 

3 

4 

5 

1 


COMMON/WORK  1/  B<35o>  .X(3^0>  iY<350»  ,YTEMP<35W)  ,A150()0>  •£  <5000  ) ♦ 
IA(5oOO>iIE(50oO)t(.A(l3o2)»i>£(2oo?)  * ICn*MI  i 3q2*2> , kjnBAS  ( j3oZ)  • 
JH(3*S0)  .ISTyPE(350)  *NAMEj20)',MTEMP(20T,CMAN,cONO,ERHAX,lFFr2. 
lNVFR0,l08Jf lRO"P»I lCHtITCHA*lTCNT,lTRFftQ*TVlN»ivOUT,JCOLP,K1NP* 
XSTAT*NRQIY*NcOL«NELfcM,NETA*NLEl.EM,NLETA»NPFLER*NGETA*NUELFM, 

NUtTA,SUMlNF,K3 

rOMMON/BLOCK/  zTOLZE»ATOLPV,ZTCoSTtNRMAX,NTMAX»NEHAX.eRO.«MA,QBA, 
QfI*GE0*08L»QPL  *OMI *UA»Q8»UC»0E*0f.Q6*0H»wj ,QL«Um»&N' QO*OR*QUtOZ 


nE  a X(IROWP) 
nP  a DE/Y«IR0WP> 

X<IR0WP)  • OP 
00  1000  I a 1*nRc“ 

IF  (I  .EQ.  IRO*P>  00  U>  1000 
RE  . X(I) 

X<I>  » DE  - Y<I>*OP 
1000  CONTINUE 
return 
pNo 


021 
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SUBROUTINE  WRETA 
C 

COMHON/MORkI/  B >350) »A O»0> tY (3s0> .VTEMP (35v) .A <5000 ) »E (5000) t 
I IA(5oOO> , IE  (5000)  »i.A  (1^02)  ,lE  <200?)  » ICN*M(1  302»*>  .KInBAS ( 1302)  . 

? JH(350) tlSTYPE(350).NAHE(20);NTEMP<20),CH4N,CONO,ERNAX.lfrEZf 

3 INVFR8»l0BJ*I«0«Pf lTcM»ITCHA»i!cNT.lTRFRO«iVIN»IVOUT*JCOtP.KINP» 
A XST  AT  »NRo*<*NCOL  tNELEN.NET  A. NLELEMfNLElA»N*»FLEM»N3ETA»NUELEM. 

5 NUETA,SUMINf$K3 

rOMMON/BLOcK/  ZlOLZEtZTULPVtZTCoSTtNRNAXtNTMAXtNEMAXtQROtBNAtOBAt 
1 QFT »CE0»QBL»GPL*QMI»8A»QB»QC.qE.QF.QQ*QH»wt»01»Qm»QN.Q0»QR,QU.QZ 

C 

NElEH  . NElEM  ♦ i 
TE(NELEM)  a IROXP 
F<NEi EM)  ■ Y(IROWP) 

C 

nO  1000  I ■ I.nHOK 
jF  (I  ,EQ.  IROttP)  00  TO  1000 
TF  t ABS(Y(I>)  .LB.  zT°LZ£)  go  t°  l0"0 
nELEM  a NElEM  «~1 
fETNELEM)  ■ I 
E*NFjEM>  * Y<I) 
looo  CONTINUE 
C 

NET*  ■ NeTa  ♦ 1 
LE(NETa*l>  ■ NfLf M ♦ 1 

return 

pNO 

SUBROUTINE  ITErOP(IPAK) 

c 

cOmmOn/MORkI/  8(330)  .A(3i>0)  tYOsO)  tYTEMP(3Su)  ,A(5000)  tE(SOOO)  t 

1 IA(5oOO> . IE (5000) .LA(li02> .LE (2002»  * 1?N*M » , 302 .2) .KIn0A$ < i 3o2> , 

2 jh (350) t iSTYpE (350) »NaME(20> tNTEMP f20) .CMiN.COND.ERMAX. IFFFZ. 

3 INVFAQ.lOBj.lRo*P.lTCH.lTCHAiiTCNT.lT$FRQMVlNtivOUT.jCOWP.KlNP. 

a xstat, nRop»ncOl. nelem, neTa,nlelem,nleTa.n*eLem,noeta,nu|lem, 

s NUETa.SUmINF.K3 

cOmmon/blOck/  ztolze » ztolpv » ztcqSt .nbmax t ntmax.nehax. qro* OMa . qba . 

1 OFI.OCO.OBL.OPL.OMI.OA.OB.QCioE.QF.OG.QH.wt «QL*Qm»QN*QO*QR,QU»QZ 

C 

IF  (IPAR  »EO.0>  GO  TO  1000 
OBJ  M-XdOBJ) 

IF  ( ?FFEZ  .EQ.  0)  OBJ  ■ SUMINF 
C 

IF (K3.NE.1)RETuHn 

WRITE (6.8000)  I TcNT.XSlAT.OBJ.lv  IN. I VOUt.CMiN, 

InETA.NELEM.TIMER 

8000  fORMAT(1H  #lS,4AA*,2X,tlP.8,4X.l6#4x.l6.AX,Ft6.B.Ax.I0,I8. 
lFfl.Z  ) 

60  TO  9000 

1000  fF(K^.NE.l)RETuHN 
WRITE (6.8100) 

BIOO  fORmaT (//8hOITcOUNT.2X6kSTATUS»4X9hO«J  VALUt.8X.5HVECIN»5X*HVEC0UT 

l.llX,2M0J.l2X,4HNETA,3Xt5HNELEM,4X,4MTIME  ) 

9000  RETURN 

end 
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o r>  o o on 


SUBROUTINE  INVfBT 

cOmmON/WORkI/  0(330)  *xoao)»r<350)  (VreM^OSut  «A(3000)  »E(50Qo>  • 

1 I*<5ooO> , IE  (3000)  .I.AU302)  »l£  (2002 t . ICn*M ( 1 3o2.2>  »KIN0A|(13O2) » 

2 Jh<350>  .ISTybEOSO)  *RAM£(2oj  ;n!emp<20)  ,$RiNtC0N6tERMAX,IrrEZ, 

3 lNVrBQ#IO0JtrRC»P*lTCH»ITcwx#iTcNT.lTRrfta«tVlN»lVOUT.JCOLR.KTNP# 
* XBT*TtNROP»NCOL»N£LEM,'YCTA,NLELEM,NLeTA#N»BLEP»R«ETAtNU£LEM» 

5 NUET  a t SUmINF  t*3 

r0MM0N/8L0CK/  Z.roL2E>^TOl.PV*2TC3Sr»NaMAA,NTnAX«NEMAX*QBO*9MA»QBA« 

1 Qr  f •«E0»Q8L*QPL»QMI«UA»Qa»QC»ije»0r  ,0G*QM»mi  »QL*aM«8N,Q0*0B,QUtQZ 

TNTeOeo  M«EG,rtPE3,VHES 

nlMfNSTON  ««EG(350)  *HHEG(33o>  *YREGti5o> 

SET  PARAMETERS 

n£Ta  • 0 
nLETA  • 0 
nGETA  • o 
wueta  • o 
mELEP  • 0 
NLELE*  ■ 0 
nGel£m  ■ 0 
hjUelEM  » 0 

nabove  ■ ® 
i.EO»  • l 
t Rl  ■ l 

*R1  • o 

Lfl*  ■ kRom  ♦ l 
««4  ■ M*OW 

PUT  SLACKS  AND  ARTIFICIALS  in  Part  4 ANu  REST  JN  PaRT  1 

nO  100  I ■ 1*NROa 

TP  ( JM ( I)  .GT.  NROP)  GO  TO  So 

I «4  ■ L R4  - 1 

mREG(LP4)  a JH  ( I ) 

VREG<LR*>  » JH ( I ) 

GO  TU  RO 

50  kRi  ■ kRi  ♦ 1 

vREG(KPl)  > JH  ( 1 ) 

9°  wREG  ( I > ■ -1 

jH ( I ) a 0 
100  CONTINUE 

kR3  > lRA  -1 
LR3  a I.RA 

nO  ?6o  I ■ LR4.AR4 
fR  ■ mbEG(I> 
mREG(Ib>»  0 

JH(IR)  ■ IR 
kIN8»SUR)  ■ IP 
200  rONT fNuE 


C 

C 


PULL  OUT  VECTORS  BeLOW  bUMP  ANu  gET  Rq»  COUNjS 


MIS  PA®  IS  BEST  QUALITY  PRACTICABLE 

jBOM  C**rY  FUKWSttSL  TO  W>,Q 


r » 


fclBNONZ  ■ KR*  - UR4  ♦ I 
iF  <K*i  ,EQ,  0)  60  TO  il90 
J ■ LR\ 

210  iV  ■ VREG ( J) 

LL  ■ L4(IV) 
kK  m L*  <IV*1)  .1 
iRcnt  • 0 
no  2f0  I ■ ll»kk 

NgNONZ  * NBNONZ  ♦ 1 
!R  • I*(I> 

TF  (HREG(Ir)  .qE.  o>  00  TO  220 
fRCNT  • IRCNT  ♦ 1 
hrE6 ( IR > • HREflllR)  - 1 
rRp  ■ jR 
229  cONTTNUE 

jF ( IRCNT  - 1)  23o»250*i00 
230  CONTINUE 

iF(K3.e0.1)PrIn!  8000 

8000  rORM«T(16HoMATRlX  SINGULAR  ) 

KlNBAs(IV)  « 0 
VREfl(J)  » VREG(Kpl) 
kR1  ■ kRI  - I 
fF  ('}  ,flT • KRI)  GO  TO  ?10 
GO  T n 210 
C 

250  vRE0<J>  ■ VREQ«KR1> 
kr1  ■ KRl  • 1 
lR3  ■ lR3  • 1 
vREQ(LR3)  ■ IV 
mRE6<LR3>  ■ I«P 
wrEG ( IBP)  ■ 0 
jM(iRP)  ■ iV 
KlNBiS (JV)  ■ JRP 
jF  ( j .AT.  NR1)  GO  TO  I1® 

60  T r>  210 

300  fF  (J  ,GE»  KRl)  GO  TO  310 
J ■ J*’ 

GO  TO  210 
C 

c pull  out  Remaining  vectors  ago*?  and  below  the 

C BUMP  AND  ESTABLISH  MERIT  COUNT*  OF  CPlUmNS 

310  nVRem  ■ 0 

TFlKRl  .EO.  0)  GO  TO  H9O 
. J ■ i R T 
32?  jv  a VPEG* 

I.L  a LA(IV) 

kK  ■ L4 ( IV*1)  . 1 

IRCNT  ■ 0 

no  8 if 0 i ■ ll»kk 

iR  ■ I* (i) 

TF(hREG(!R)  ,Ne,  -2)  GO  TO  *00 
C 

C PIVOT  above  BUMP  IPaRT  OF  L) 

C " 

NABOVE  ■ NABOVp  * 1 
iROWo  ■ IR 
CALL  UNPACK (IV) 


oo  O O oo  oo  o 


t 


SHIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
XQOM  CO t*  FUKtUSHB)  TO  DDC  


call  bret* 

nLETA  ■ NET* 

JH(IR)  a IV 
kINBAS(IV)  ■ IR 
vR£6<J>  • VREG(KRl) 
kR1  ■ KRl  • * 
nVReM  ■ NVREH  ♦ 1 
hRE<3(IR>  • IV 
60  TO  9*0 
C 

*00  tF  (HREQ(Xr)  .Gfc,  o)  60  JO  800 
fRCNf  a IRCNT  ♦'! 

A fRP  > IR 
80?  cONTTNUE 

c 

IF  (IRcNT  - 1)  810*900*1000 
81 0 CONTINUE 

tF(K3.E°.1)RrIn!  8000 
KINBaS(IV)  • 0 
VREG«J>  ■ VREG(KPl> 
nVRem  a NVR£M  ♦ 1 
*R1  ■ kRI  • 1 

IF  ( »OT*  KRI)  GO  TO  ioio 
r,0  Tft  320 


Pu!  VECTOR  BELO*  BUMP 

900  VR£G(J>  ■ VREG(KRl) 
nVREm  ■ NVREH  ♦ 1 
kRI  ■ KRl  • 1 
lR3*  lR3  • 1 
VREG(LR3)  ■ IV 
mREG<L«3)  ■ IHp 
hREO(IHP)  s o 
jH(iRp)  a IV 
KlNBAS(IV)  ■ IRP 


change  pom  counts 

9*0  f)0  950  II  ■ LL.KK 
I IR  b I*  ill) 

tF  (HREG(IiR)  ,6E « 0)  «°  TO  950 
hREGMIR>  • HReGiIIR)  * i 
950  CONTINUE 

IF  (J  ,GT.  KRI)  fcO  TO  lOlO 
. GO  TO  32« 

loo?  tF  (J  .GE*  kRI)  GO  TO  1010 

.)  ■ J*’ 

GO  TO  120 

LOlO  IF  (NVREH  .GT.  0)  GO  JO  310 

geT  merit  counts 

1020  tF  (KR)  ,EQ.  0)  60  TO  1190 

oo  1T00  j b lRi»kRi 

IV  a VREG  < j) 
lL  . L* ( IV) 

KK  b LA ( IV* 1 ) • 1 
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IBIS  PAGE  IS  BEST  QUAI/ITT  PRACTICABLE 
fPfM  COfc'  Y FURBISHED  TO  SDR 


tmcnt  • 0 

t)0  1*50  I ■ LL«Kk 

lR  ■ I A < I > 

IF  (HRfG ( IR)  ,g£,  0)  60  TO  1050 
TMCNT  » IMcNT  - (HREGHRJ  ♦ !> 
*050  continue 

hREG(J)  ■ IMCNT 
1 l00  cONTtNUE 


c 

c 

c 

c 


, . tSO  • 1 

H06  <*Rl  ,LT, 

ISO  ■ **ISD 
fiO  TO  1106 
1106  fSO  • ISO  - 1 


Sort  colomns  into 
using  shell  sort 

2*150)  GO  TO  1108 


merit  oroe1? 


c ENO  OF  INITIALIZATtON 

llol  IF  <ISo  .Lg.  0)  GO  TO  HOT 
, Is*  • l 
no?  tsj  ■ isk 

iSL  * ISK  ♦ ISO 
tSY  ■ MREG(ISL) 

, . fsz  » VREG (ISL) 

lloj  jF  <fsv  #LT»  MfiEG(ISJ))  60  TO  ll04 
HOS  iSL  . ISJ  ♦ ISO 


mREG<ISL>  ■ 1ST 
VREG  (ISL)  ■ ISZ 
iSK  ■ ISK  ♦ 1 

TF  (<ISK  ♦ ISO)  ,LE.  KMl)  GO  TO  1102 
(SO  > (ISO  - 1)  / 2 


GO  T«  HOl 

1104  TSl  ■ ISJ  ♦ ISO 

mReq(ISL)  » MRtG(ISj) 
VREG<ISL)  ■ VRpG(ISj) 
iSj  ■ ISJ  - ISd' 

IF  (ISJ  .GT*  0)  gO  TO  1103 
C,0  To  1105 
HOT  CONTINUE 


ENO  OF  SONT  routine 


PUJ  OUT  BELUN  BUMP  ETAS  (PART  Up  U> 


1190  nSlcK  ■ o 
nBELOM  m 0 

nelast  ■ nemax 
nTLAST  ■ ntmax 
I.E(NTLAST  ♦ 1) 


■ NELAST  ♦ 1 


LH  ■ LR3 

IF  (LRT  *Ge*  LR4)  lR  ■ L*4 
IF  (LR  .OT.  KR*)  GO  TO  «OSO 
jK  a KR4  ♦ 1 

00  2000  JJaLR t kR* 
jK  ■ JK  - 1 

IV  . VREG(JK) 

1 • MREG(JK) 
nOElOM  • N8EL0W  * 1 


1 


t 


i 
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r»or»  o o o r> 


t 


THIS  PAGE  IS  BEST  QUALITY  FRACTICABI* 
from  copy  furnished  ioddc 


IF  (IV  .6T,  NHO")  00  TO  1200 
nSlc*  ■ nSlc*  ♦ 1 

1200  Lc  . |_ A ( I V > 

KK  ■ LA ( X V* 1 ) -i 
IF  (KK  ,QT.  LL)  GO  TO  *?®0 
125?  TF  ( «B$  ( A <LL>  - 1.)  . Lfc.  ZTOLZE)  GO  TO  ?0?0 

1300  kjUETA  • NUETA  ♦ 1 
r>0  1*00  J a LL*KK 
T«  ■ I A ( J 1 

IF  ( IR  .eO.  I)  00  TO  UoO 
TE<NrLAST>aIR 
F<NELAST)  ■ A ( j) 
nEL AST  ■ NELAST  - 1 
NUELEH  ■ NuELEw  ♦ 1 
GO  Tft  U00 
1390  rF  ■ A ( J) 

1A00  rONTINllE 

lE  (NfLaST)  ■ 1 
f (NEi AST)  ■ EP 
l E (nTL AST)  ■ NELAST 
NELAST  • NELAST  • 1 
nTiaST  a NTLAST  • 1 
nUelEM  ■ NUELEH  ♦ 1 
2000  ^ONTINLE 

2050  tF(KP1  .E<3.  0>  00  TO  350? 

00  L-U  DECOMPOSITION  Of  B^MP 

nO  300"  J ■ L«l»KRl 
jv  • vosg < j) 

CALL  UkPACK(IV) 

CALL  FTRAN(2» 

tRowp  ■ o 

IRCmtN  ■ "999999 
DO  2*00  I a 1»N«0« 

rF  ( APS ( V ( I ) ) .LE.  ZTOLPV)  GO  TO  21 "0 
rF  (HBEG(I)  .GE.c)  go  lo  2100 
rF  (MRFGU)  #Le»  IRcmIN)  GO  TO  jlOO 
I«CMIN  • MREG<l> 
jROWp  * I 
2100  CONTINUE 

IF  (IROWP  .GT.  o)  GO  TO  215o 
ifikT.eo. Sprint  aooo 
kInbIs ( IV)  ■ 0 

GO  TO  3000 

2150  fNCR  ■ HREgURO«P)  * 3 

WRITE  L ANO  U ETAS 

TF  (J  .EO*  KRl)  gO  TO  2160 
NELEN  ■ NELEN  ♦ I 
iE(NrLEM)  • IRQWP 
F<NElEM)  ■ T(IROhP) 

2160  no  2t00  i ■ 1»nRo* 

• IF  ( j .EQ*  IROmP)  GO  TO  2300 

rF(  ABS(V(D)  ,LE»  ZTOlZE)  GO  To  23(1'* 


1 


■ 

i 


f 


Cr-2J 


«IS  W®  IS  O'”*1”’ ™“"T“ 

J5J0H  OOlFY  TUBWlSHiiD  10  UOfi  — ^ 


iP  (HRpOd)  tOE*  0)  GO  TO  2200 
C 

C L ETA  ELEMENTS 

C ... 

nELEM  • NELEM  « ] 
fE (NrLfM)  • I 

rlNEl  EM)  ■ Yd) 
pO  TO  ?300 
C 

C U £TA  ELEMENTS 

C 

2200  fE(NELAST)  ■ I 
F<NEi  AST)  • Y(D 
n|LA«T  ■ NEL^ST  - 1 
„ wUClEM  • NuELEm  ♦ i 
?300  CONTINUE 

c 

jH(IROnP)  ■ IV 
kINbAS(IV)  ■ ISOwP 
nUeta  ■ NUgTA  ♦ t 
tE(NfLaST)  ■ IROvP 
TP  ( 1 ,NE»  KRl)  GO  TO  £330 
F<NELAST)  ■ Y(lNo»P) 

GO  TO  2340 
2330  F(NELAST)  ■ l. 
nETA  ■ NETA  ♦ l 
lE(NETA#I)  ■ NELEM  ♦ 1 
?340  NUELEM  • NuelEM  ♦ I 
LE(NTLAST)  ■ NeLaST 
nELAST  • NELAST  - i 
A)TlaST  • NTLAST  - 1 
C 

C UPOaTe  RON  COUNTS 

C 

oO  235"  I ■ 1»N«0* 

fP  < A0S(Y(I))  «lE.  ZTOlZE)  GO  TO  23*0 
IP  (HREG(Ij  ,GE.  0)  GO  TO  2350 
hReg‘I>  ■ HREG(I)  • INCR 
I P (HREGd)  *GE«  °)  MKEGd)  » *1 
2350  CONTINUE 

H«EG(IROnP)  a o 
3000  CONTINUE 
C 

c merge  l and  u etas 

c ' 

3500  NLETA  ■ NETA 

nETA  ■ NLETA  ♦ NUETA 
nLELEM  ■ NELEM 
NELEM  ■ NLELEM  ♦ NUELEM 
TP  (NUELEM  .E0*  0)  GO  TO  3550 
C*LL  SHPTE 
C 

C INSERT  SLACKS  POP  DELETED  COLONS 

C 

3550  nO  3600  I ■ 1»nRC* 

IP  (JHd)  ,NE.  0)  GO  TO  3*00 
JH(I)  • I 

tnonp  ■ I 


non  o <"»  ooo  ooo 


IHIS  PAGE  IS  BBS!  QUALITY  FRACXICASLS 

roan  copy  ju^lshjEd  ro  dd,g  


CALL  UNPACK!  I) 
CALL  FTRAN(l) 

call  breta 
3600  CONTINUE 


UPUaTe  X 

rALL  SnIFTr ( 1*3) 

CALL  FTRAN(l> 

CALL  ShIFTR<3*2) 

PHInT  STATISTICS 

nOFO  ■ NELEM  - NeTa 
mSTP  ■ NRON  - nSlCK 
tF(k3.nE.1)HETuRN 

W^ITf (6*500)N8NON2tNSTR,NABOVE»NBELOw*NLELE«.NLETAtNUELEM*NUETA» 
1NOF0,NETA 

SoO  pORmAT(18HoINVeRT  STATISIICS/1H  ,U,i4H  NONA  IN  BAsiS/lH  »U, 

1 ?8H  STRUCTURAL  COLUMNS  IN  BASIS/1M  *TA»19H  VICTORS  ABOVE  BUMP/lH  * 
?T*»19H  VECTORS  BELOB  »UmP/3H  LI,I5»5h  N0NZ.i5,5H  ETAS/3M  U|,I5# 

3SH  nONZ » l5»5H  E T aS/8h  JOT ALSI » I5» Uh  OFF  OI*ft  NOn2,I5.5H  ETA5  ) 

RETURN 

fNq 

SUBROUTINE  shfte 


cOhmON/BORkI/  8<35q>  *A(3&0) *Y (350»  »YtEMP(35u> »A<5000>  »E<5000> t 
1 I A <5p00>  » IE (50o0) ,L* (1^02) *l£ (200?) * ICn*M(1 302*2) *KInBA$ ( 130?) , 

? JH(350) * ISTtPE ( 350 ) *NAME<20) »nTEMP  *20) ,CMlNtCCNO,£RMAX* JFFeZ* 

5 INVFRO*IOBJ»lHO*P*I «CH»ITCHA»lTCNT.lTBFBO»TVlN*ivOUT*JCOLp,KINPt 
A XSTAT»NRoB*NCOL»NELEMtNETA.NLELEM,NLCTA»N»frLEM*NOETA,NUELEM* 

S NUETA*SUMINF,K3 

COMMON/BLOCK/  ZToLZE*Z10LPV.ZTCoST*NOMAX,NTh*x»NEMaX.QRO*«MA*QBA* 

1 OF!*CEO*QBL»OPL*QMI»UA*Qb*UC*oE*OF.OG»QH*WT*OLfQM*QN«C)0«QR,OUfOZ 

Sm.IfT  IE  AND  E Of  g ele«en!s 


* 


♦ 


» 


» 


nF  ■ NEMAX  - NuELEM  ♦ 1 
|NCR  * 0 

no  1000  1 . nf.nemax 
jncr  ■ inCr  ♦ 1 
lE(NLELEM  ♦ INCH)  ■ IE ( I ) 

E (NLfLFM  ♦ INCR)  ■ Ell) 

1000  CONTINUE 
C 

toif  • nemax  - nlelem  - nuelem 

nF  a NTMAX  - NUETA  ♦ I 
iNcR  ■ 0 

pO  2000  I ■ NF. NTMAX 
INCR  ■ INCR  ♦ 1 

lE(NLETA  ♦ INCR)  a LE ( I ) - IOIF 
2000  CONTINUE 

I E<NETAM)  ■ NeLeM  ♦ 1 

RETURN 

rNO 

SUBROUTINE  GETPHI ( I * J* A.F ) 

RETURN 

fno  C-29 


